PRIMERA CLASE: REGRESIÓN LINEAL
Maestría en Investigación Operativa y Estadística
INTRODUCCIÓN
DOCENTES
Héctor Hernán Montes García.
Ingeniero industrial de la Universidad Tecnológica de Pereira.
Estudios en Maestría en Ciencias con orientación en Matemáticas
Desempeño laboral
Científico de datos con más de 3 años de experiencia en la construcción de modelos de aprendizaje automático y soluciones de minería de datos para obtener mejores conocimientos comerciales. Experiencia utilizando SQL, bibliotecas de Python (matplotlib, seaborn, flask, scikit-learn, pandas, numpy), modelos matemáticos y estadísticos para ofrecer soluciones robustas que agregan valor al negocio.
Ha trabajado para clientes del sector industrial y financiero en México y Colombia, en áreas tales como:
- Mantenimiento predictivo
- Diseño de campañas comerciales basadas en datos
- Reconocimiento de imágenes
- Modelos de procesamiento de lenguaje natural (NER).
Julián Piedrahíta Monroy
Ingeniero industrial.
Magíster en desarrollo agroindustrial
Universidad Tecnológica de Pereira
Desempeño laboral
- Consultor en soluciones análiticas para instituciones educativas.
- Analista de datos en el observatorio social de la UTP.
- Actualmente diseñador de tableros informativos (dashboards) para la Universidad Pedagógica Nacional de Bogotá y la misma Universidad Tecnológica de Pereira.
- Docente catedrático en el área de informática y estadística general.
- Uso principal del lenguaje R y sus librerias Tidyverse, RMarkdown, RShiny y otras.
Contenido del curso
- T1. Regresión Lineal, suposiciones y requisitos.
- T2. Estimación de los parámetros e interpretación de los valores.
- T3. Pruebas de hipótesis relacionadas con los parámetros.
- T4. Evaluación de modelos de acuerdo a las suposiciones de normalidad e independencia.
- T5. Transformaciones de las variables, y sus implicaciones en las pruebas de hipótesis.
- T6. Series de Tiempo. Estacionariedad y cómo obtener una serie estacionaria a partir de una que no lo es. Correlogramas.
- T7. Modelos de Box y Cox. Estimación de parámetros.
- T8. Criterios de evaluación de un modelo de series de tiempo. Transformaciones.
Bibliografía.
Regression Modeling and Data Analysis with Applications in R. Samprit
Chatterjee, Jeffrey S. Simonoff
Montgomery, D. C., Peck, E. A., & Vining, G. G. (2002). Introducción
al análisis de regresión lineal (3ª ed.). Limusa Wiley.
Time Series Analysis and Its Applications: With R Examples by Shumway
and Stoffer.
Time Series Analysis: Forecasting and Control by Author(s): George E. P.
Box, Gwilym M. Jenkins, Gregory C. Reinsel, Greta M. Ljung.
https://fhernanb.github.io/libro_regresion/rls.html#modelo-estad%C3%ADstico
https://bookdown.org/victor_morales/SeriesdeTiempo/
Motivación
Los métodos de regresión son técnicas estadísticas utilizadas para modelar y comprender la relación entre variables. La motivación detrás de estos métodos radica en la necesidad de analizar y cuantificar cómo una o más variables predictoras influyen en una variable respuesta. Los métodos de regresión son ampliamente utilizados en diversos campos, incluyendo la investigación científica, la economía, la medicina, la ingeniería y más. A continuación, se detallan algunas de las principales motivaciones detrás de estos métodos:
Modelado de Relaciones:
La motivación fundamental de la regresión es modelar la relación entre variables. En muchos casos, existe la necesidad de entender cómo una variable dependiente cambia en función de una o más variables independientes. Por ejemplo, en la economía, podría ser necesario entender cómo las tasas de interés afectan el gasto del consumidor.
Predicción y Estimación: Los modelos de regresión también se utilizan para hacer predicciones y estimaciones. Dado un conjunto de datos históricos, los modelos de regresión pueden ayudar a predecir valores futuros de la variable dependiente en función de los valores de las variables independientes. Esto es especialmente útil en campos como la meteorología, la finanzas y el marketing.
Control y Optimización: En algunos casos, se utilizan modelos de regresión para optimizar y controlar procesos. Por ejemplo, en la manufactura, los modelos de regresión pueden usarse para identificar las condiciones ideales que conducen a la máxima eficiencia o calidad del producto.
Validación de Teorías: En la investigación científica y social, los modelos de regresión pueden usarse para validar o refutar teorías existentes. Al comparar los resultados del modelo con las expectativas teóricas, se puede evaluar la validez de las hipótesis.
Gestión de Riesgos: Los modelos de regresión también se utilizan para evaluar riesgos y tomar decisiones informadas. En finanzas, por ejemplo, se pueden utilizar para evaluar el riesgo de inversión y la exposición a diferentes factores del mercado.
Ejemplo
El modelo CAPM (Capital Asset Pricing Model) es un modelo de valoración de activos financieros desarrollado por William Sharpe que permite estimar su rentabilidad esperada en función del riesgo sistemático. Su desarrollo está basado en diversas formulaciones de Harry Markowitz sobre la diversificación y la teoría moderna de portafolios. Dos puntos claves son:
- El modelo CAPM es utilizado para calcular la rentabilidad que un inversionista debe exigir al realizar una inversión en un activo financiero en función del riesgo que está asumiendo.
- El modelo CAPM establece una relación lineal entre el rendimiento esperado de un activo y su riesgo sistemático, medido por su \(\beta_i\).
La fórmula del modelo CAPM es la siguiente:
\(E(r_i)=r_f+\beta_i*(E(r_m)-r_f)\)
Donde:
\(E(r_i)\) es la tasa de rentabilidad esperada de un activo concreto. \(rf\) es la rentabilidad del activo sin riesgo. \(\beta_i\) es la medida de la sensibilidad del activo respecto a su benchmark. \(E(r_m)\) es la tasa de rentabilidad esperada del mercado en que cotiza el activo.
Análisis de Causa y Efecto: Los modelos de regresión pueden ayudar a establecer relaciones causa-efecto entre variables. Esto es útil para comprender cómo los cambios en una variable influyen en otras variables y viceversa.
En resumen, la motivación detrás de los métodos de regresión radica en la necesidad de entender, modelar, predecir y cuantificar las relaciones entre variables en una variedad de campos. Estos métodos permiten tomar decisiones informadas, realizar análisis profundos y obtener información valiosa a partir de los datos disponibles.
1. INTRODUCCIÓN NO FORMAL AL MODELO LINEAL SIMPLE
En este capítulo introduciremos el modelo lineal simple partiendo de datos reales del precio del aguacate y el precio del dólar. Más que una presentación formal o matemática del modelo, nuestro objetivo será construirlo paso a paso usando las ideas estadísticas detrás del método. Confiamos en que esto le al lector una idea más precisa de cómo éste modelo usa los datos dados para relacionar las variables de inteŕes, y que comprenda sus fortalezas pero también sus limitaciones.
1.1 Lectura de datos
A continuación, cargaremos dos tablas de datos que contienen información sobre el precio del dolar, también conocida como tasa representativa del mercado, y otra con la información del precio del aguacate Hass tomada de una fuente gubernamental y de la página del SIPSA.
Los enlaces para acceder a información relacionada son:
https://www.agronet.gov.co/estadistica/Paginas/home.aspx?cod=11 https://www.bde.es/webbe/es/estadisticas/temas/tipos-cambio.html
1.2 Creación de funciones.
Se programa una función de flextable que mejorará la impresión de tablas de resumen. Para usarla sólo bastará invocar la función ftable() al final de la sentencia (usando %>% tidyverse) o con el ftable() encapsular o encerrar la tabla que se quiere ajustar.
# Funciones
## Función para crear flextable
ftable <- function(x) {
x %>%
flextable() %>%
theme_vanilla() %>%
color(part = "footer", color = "#666666") %>%
color( part = "header", color = "#FFFFFF") %>%
bg( part = "header", bg = "#2c7fb8") %>%
fontsize(size = 11) %>%
font(fontname = 'Calibri') %>%
# Ajustes de ancho y tipo de alineación de las columnas
set_table_properties(layout = "autofit") %>%
# width(j=1, width = 3) %>%
align(i = NULL, j = c(2:ncol(x)), align = "right", part = "all")
}Como ejemplo del uso de la función:
Sin flextable
## valor unidad vigenciadesde vigenciahasta
## 1 4,761.64 COP 22/12/22 22/12/22
## 2 4,781.28 COP 20/12/22 20/12/22
## 3 4,802.48 COP 17/12/22 19/12/22
## 4 4,836.24 COP 13/12/22 13/12/22
## 5 4,815.99 COP 10/12/22 12/12/22
Con flextable
valor | unidad | vigenciadesde | vigenciahasta |
|---|---|---|---|
4,761.64 | COP | 22/12/22 | 22/12/22 |
4,781.28 | COP | 20/12/22 | 20/12/22 |
4,802.48 | COP | 17/12/22 | 19/12/22 |
4,836.24 | COP | 13/12/22 | 13/12/22 |
4,815.99 | COP | 10/12/22 | 12/12/22 |
1.3 Depuración de los datos
Dolar
Como en todo proceso de análisis de datos, es necesario realizar una depuración y ajuste de los datos. Para este caso, fue necesario trabajar sobre la variable de la fecha y el valor.
# Vamos a sacar un valor promedio de cada variable por mes
#Para el dolar vamos a tomar el campo vigenciadesde
dolar <- dolar %>%
mutate(fecha = as.Date(vigenciadesde, format = "%d/%m/%y")) %>%
mutate(mes = month(fecha), anio = year(fecha)) %>%
mutate(valor = as.double(str_replace(valor,",",""))) %>%
group_by(anio,mes) %>%
summarise(precio_dolar = mean(valor), .groups = "drop") A continuación se muestra una fracción de los datos depurados.
Aguacate
mercado | producto | fecha | precio_kg |
|---|---|---|---|
Bogotá, D.C., Corabastos | Aguacate hass (Semanal | Mensual) | Sabado, Junio 23 de 2012 | 3,483 |
Bogotá, D.C., Corabastos | Aguacate hass (Semanal | Mensual) | Sabado, Junio 30 de 2012 | 3,459 |
Bogotá, D.C., Corabastos | Aguacate hass (Semanal | Mensual) | Sabado, Julio 7 de 2012 | 3,478 |
Bogotá, D.C., Corabastos | Aguacate hass (Semanal | Mensual) | Sabado, Julio 14 de 2012 | 3,493 |
Bogotá, D.C., Corabastos | Aguacate hass (Semanal | Mensual) | Sabado, Julio 21 de 2012 | 3,470 |
# 🖇️ Se genera el vector de meses para usarlo más abajo con la función match.
meses <- c("Enero", "Febrero", "Marzo", "Abril", "Mayo",
"Junio", "Julio", "Agosto", "Septiembre", "Octubre",
"Noviembre", "Diciembre")
hass <- hass %>%
# Otra forma de sacar el anio.
#mutate(anio = substring(fecha,nchar(fecha)-4,nchar(fecha))) %>%
#mutate(espacio = grepl(", ",fecha))
mutate(mes = sapply(strsplit(fecha, " "), "[", 2),
anio = sapply(strsplit(fecha, " "), "[", 5)) %>%
mutate(anio = as.double(anio)) %>%
# 🖇 Se utiliza la función match con el vector meses.
mutate(mes = match(mes,meses)) %>%
group_by(anio,mes) %>%
summarise(precio_aguacate_kg = mean(precio_kg), .groups = "drop") 1.4 Diagrama de dispersión.
Cómo se puede notar uno intuye que hay alguna relación entre ambas variables en la medida en que parece suceder que cuando el precio del dólar aumenta, también lo hace el precio de aguacate. La relación no es perfecta pues no vemos una curva suave que sea capaz de pasar por todos los puntos, más bien debemos reconocer que hay ciertas variaciones en el proceso. Estas variaciones pueden deberse a factores no contemplados por el modelo, después de todo no podemos esperar que el precio del aguacate dependa exclusivamente del precio del dólar.
1.5 Generando un modelo súper simplificado (modelo ingenuo)
Ahora vamos a generar el modelo más sencillo que podamos imaginar para predecir el precio del aguacate, al cual llamaremos un modelo ingenuo. La razón del apelativo es que el modelo “ingenuamente” supondrá que es posible predecir el precio que tendrá el aguacate con el promedio histórico.
La siguiente gráfica muestra la linea horizontal que simboliza el promedio histórico del precio del aguacate. También debe notar que predecir basado en esta línea es despreciar cualquier aporte de la variable precio del dólar en la predicción, es decir, es un modelo que no explota la relación (sea esta débil o sea esta fuerte) entre el precio del aguacate y su variable explicativa precio del dólar.
hass_dolar %>%
ggplot(aes(x= precio_dolar, y= precio_aguacate_kg )) +
geom_point()+
geom_hline(yintercept = mean(hass_dolar$precio_aguacate_kg), color = "red")+
geom_text(aes(x= 5700, y = mean(hass_dolar$precio_aguacate_kg),
label = paste("Precio promedio histórico: $",
round(mean(hass_dolar$precio_aguacate_kg), 2))),
hjust = 1.2, vjust = -0.2, color = "red"
) +
theme_light()Observe que el modelo ingenuo no sería tan inadecuado en dos circunstancias:
Si los datos se encuentran muy cercanos al precio histórico de manera consistente sin importar en que valor del precio del dólar me posicione. Lo que vendría a indicar un escenario donde moverme a lo largo de diferentes valores del precio del dólar no genera ningún patrón de distanciamiento frente al valor histórico de referencia.
Si los datos del precio del aguacate se alejan de la línea histórica de referencia pero el precio del dólar tampoco puede seguirlos, es decir, si la relación entre precio del dólar y precio del aguacate se parece a una nube de puntos sin ningún patrón o tendencia evidente. En este caso cualquier cosa que supongamos de la relación entre precio de dólar y precio del aguacate será producto de nuestra imaginación 🤦🤦
En estos dos escenarios es mejor retener el modelo ingenuo a falta de una variable que de verdad explique lo que pasa con el precio del aguacate.
Ahora bien, ninguno de los dos escenarios anteriores parece ser el nuestro, más bien acá se nota que el precio del dólar si está acompasado con el precio del aguacate. Sin embargo antes de entrar en la búsqueda de esas relaciones es importante notar dos cosas:
El modelo ingenuo es mi modelo de referencia, esto significa que si un modelo predictivo construido para predecir el precio del aguacate es mejor, lo tendrá que ser respecto a este modelo ingenuo.
Es posible medir el desajuste actual de mi modelo ingenuo tomando la distancia de cada punto a la recta horizontal de promedio histórico, esto es:
1.6 Calidad de ajuste del modelo ingenuo
Sea \(e_i = y_i - \bar{y}\) las diferencia entre el dato del precio de aguacate \(y_i\) y el promedio histórico \(\bar{y}\). Como tengo muchas diferencias, será necesario definir una métrica resumen, por conveniencia elijamos la suma de los cuadrados de las diferencias, debido a que si sumamos las diferencias, éstas se me compensarán, pues diferencias negativas cancelarán las positivas, y no quiero esto. Más bien me interesa que ambas sumen a mi medida de desajuste. Llamemos a esta medida la suma de cuadrados del error:
\[SCE_{\text{ingenuo}} = \sum^{n}_{i=1}e_i²=\sum^{n}_{i=1}(y_i - \bar{y})²\]
En la definición de la cantidad quisimos usar el subíndice “ingenuo” para nombrar al SCE, esto con el fin de hacer énfasis que tal métrica se asocia al modelo, si el modelo cambia a otro tipo de modelo el SCE cambiará también: será la distancia de cada \(y_i\) a la curva definida por el otro modelo. Veremos esto más adelante. No obstante ya podemos sacar nuestras primeras conclusiones:
- Entre más grandes sean cada una de las distancias \(e_i\) individuales, más grande será SCE.
- Todas las distancias contribuyen con igual importancia al SCE, excepto en lo que se refiere a su magnitud no hay porque pensar que un distanciamiento de \(u\) unidades sea más importante que otro de las mismas \(u\) unidades. Esto puede ser obvio pero no lo es cuando, por ejemplo, queremos castigar más las distancias que se dan en la zona central del gráfico que las que se dan en los extremos del gráfico.
- El SCE es sensible a la cantidad de datos, razón por la cual entre más datos (es decir entre más grande sea n), más grande será SCE. Esto dificulta seriamente la posibilidad de comparar modelos que fueron calculados sobre una cantidad distinta de puntos.
De acuerdo con lo anterior modifiquemos un poco nuestra definición de desajuste, así:
\[MSE_{\text{ingenuo}} = \color{red}{\frac{1}{n}}\sum^{n}_{i=1}(y_i - \bar{y})²\]
Es decir, nuestra medida de desajuste será el promedio de la suma de las diferencias cuadradas de cada dato \(y_i\) respecto a su media \(\bar{y}\), dicho de otra manera la varianza de la variable \(Y\) a predecir!!!!
A continuación invitamos al estudiante a calcular la varianza de la variable a predecir, es decir, el desajuste del modelo ingenuo.
# En esta sección el estudiante desarrollará los cálculos.
# 💡 Recordemos que la varianza de los residuales respecto del modelo ingenuo, es la misma varianza de la variable a predecir.Observe el siguiente gráfico en el que hemos pintado las distancias que hay entre cada dato \(y_i\) observado y el promedio histórico (valor predicho \(\bar{y}\) por el modelo ingenuo)
hass_dolar %>%
ggplot(aes(x= precio_dolar, y= precio_aguacate_kg )) +
geom_point()+
geom_hline(yintercept = mean(hass_dolar$precio_aguacate_kg))+
geom_segment(aes(xend=precio_dolar, yend=mean(precio_aguacate_kg)),
col='red', lty='dashed')+
theme_light()1.7 Construyendo un modelo lineal simple con los datos
Este es el momento de comenzar a usar la información extra con la que contamos, esto es, el precio del dólar, el cual se convertirá en una variable predictora del precio del aguacate. En estadística no hay una única forma de hacer esto, pero lo usual es partir de especificaciones muy sencillas. Veamos:
Sea \(Y\) el precio del aguacate y sea \(X\) el precio del dólar. Investiguemos si este modelo sencillo nos funciona:
\[y = \beta_{0} + \beta_{1} * x + e \text{ con } e \text{~}N(0, \sigma^2)\]
1.8 Algunas observaciones sobre las fortalezas y limitaciones del modelo
Expliquemos qué implica la ecuación anterior haciendo ciertas anotaciones de gran relevancia práctica para entender cómo los estadísticos piensan al momento de plantear modelos:
Observación genial 1:
Para un x fijo, es decir un x dado, es posible calcular el valor de \(y\) con el modelo anterior, donde tácitamente estamos diciendo que el valor de \(y\) dependerá de \(x\) en forma exacta, excepto por una perturbación aleatoria \(e\) que representará la parte del valor de \(y\) que no puede ser capturada por el término \(\beta_{0} + \beta_{1} * x\).
¿Cuándo sucederá que el valor de \(y\) pueda ser exactamente capturado por el término \(\beta_{0} + \beta_{1} * x\)? Cuando \(y\) dependa en forma lineal exacta del valor de \(x\). Esta es una condición demasiado fuerte que rara vez ocurre en la práctica, por eso agregamos la perturbación \(e\) como una forma de modelar la incertidumbre en la determinación de \(y\) dado un valor de \(x\). Esto no soluciona todos los problemas pero ayuda a obtener un modelo relativamente más realista.
¿Cuáles son los riesgos que implica asumir a priori un comportamiento específico para el error?
Observación genial 2:
La perturbación \(e\) puede recibir varios nombres en estadística, la encontrarás nombrada como: perturbación, choque, error, o residual. Lo importante no es el nombre que reciba, sino las condiciones que vamos a suponer sobre sus valores.
En particular no vamos a exigir que \(e\) tenga valores predecibles, muy al contrario vamos a exigir que sea un valor aleatorio, precisamente porque está modelando incertidumbre. Pero esto no implica que no podamos poner ciertas condiciones sobre el tipo de aleatoriedad deseada para \(e\). En este caso supondremos que \(e\) se distribuye como una variable aleatoria normal con media \(\mu=0\) y desviación estándar \(\sigma\). Conviene aclarar que esto es una suposición, no implica que efectivamente los errores vayan a cumplir tal supuesto.
Observación genial 3:
¿Por qué suponemos normalidad para \(e\)? Porque si \(e\) es en efecto una especie de componente incierto en el valor de \(y\), entonces muy seguramente será la suma de muchos efectos independientes que lo están provocando. En nuestro contexto estos efectos inciertos pueden ser:
- \(F_1\): Cambios en la productividad de los cultivos.
- \(F_2\): Condiciones climatológicas.
- \(F_3\): Precios de los insumos agrícolas.
- \(F_4\): Capacidad de negociación de los productores.
- \(F_5\): Presencia o no de subsidios estatales.
- \(F_6\): Precios de productos complementarios o sustitutos.
…
- \(F_k\): precio de la gasolina.
Y un gran etcétera.
Es decir, existen innumerables factores, distintos al precio del dólar que impiden que el precio del aguacate pueda ser determinado de forma exacta sólo por observar el valor del dólar. En estadística se ha estudiado que cuando no estamos midiendo los demás factores que afectan el valor de una variable, y dejamos que estos contribuyan de forma desacoplada e independiente a dicho valor, el efecto general \(e\) que resume el efecto global de todos los factores a la vez, suele comportarse bajo la distribución normal sencillamente porque hay un teorema en matemáticas, conocido como el Teorema del Límite Central que básicamente así lo garantiza, al postular que: la suma de \(k\) efectos aleatorios independientes tiende a adoptar el comportamiento normal conforme la cantidad k de efectos crece.
Es importante destacar que el Teorema del Límite Central tiene algunas condiciones y suposiciones, como la independencia de las variables aleatorias y que la suma debe efectuarse sobre una cantidad \(k\) de ellas suficientemente grande. Éste es uno de los conceptos fundamentales en estadística y es ampliamente utilizado en la teoría y la práctica de esta disciplina.
Observación para nada genial 4 😥😢:
Excepto por el componente aleatorio \(e\), el modelo restringe la relación de \(Y\) y \(X\) al universo de relaciones lineales. Existirá una posible relación por cada par de valores \(\beta_0\) y \(\beta_1\) que decidamos elegir. Pero por más que nos esforcemos en modificarlos siempre conducirán a relaciones representadas por líneas rectas, de ahí que el modelo asuma el nombre de regresión lineal.
Si queremos capturar otras posibles dependencias de \(Y\) respecto a \(X\) podemos generalizar la relación usando \(y = f(x) + e\) donde \(f(x)\) puede ser una nueva especificación funcional con otra estructura deseada, por ejemplo un polinomio o cualquier otra función que querramos definir. Lo importante es que detrás de la definición de \(f(x)\) haya alguna justificación producto de haber analizado los datos en búsqueda de relaciones que capturen bien la dependencia.
Observación genial final:
Son los datos observados para cada valor de \(y\) y \(x\) los que nos deben informar sobre el par de valores \(\beta_0\) y \(\beta_1\) que crean la relación lineal que mejor representa nuestra nube de puntos, pero de nada sirven los valores observados si no definimos una métrica que nos informe sobre la calidad del ajuste.
1.9 Introducción a las medidas de ajuste para modelos no ingenuos
De la observación final del apartado anterior nos queda una inquietud: ¿Cómo decidimos el par de valores a asignar a los parámetros del modelo lineal si no tenemos un criterio para poder saber cuál es el mejor modelo?
Para salir de este embrollo, reconozcamos que nunca en un problema real sujeto a incertidumbre una recta pasará exactamente por todos los puntos, razón por la cual aparecerán componentes de error \(e_i\) por cada \(y_i\) y \(x_i\) observado. Vimos que esto ocurrión con el modelo ingenuo y por supuesto ocurre también para un modelo no ingenuo. La idea es entonces reducir la magnitud de estos errores al mínimo posible y el par de valores \(\beta_1\) y \(\beta_2\) que así lo logren será nuestra elección óptima de cara al objetivo de minimización del error.
Otro criterio de ajuste que se podría usar es el de maximizar la probabilidad de ocurrencia de los valores observados bajo el supuesto de que el modelo usa valores \(\beta_0\) y \(\beta_1\) previamente especificados. De esta manera si un par de valores hace menos probable haber obtenido nuestros datos observados deberán descartarse.
¿Pero cómo podemos medir la probabilidad de ocurrencia de nuestros datos observados una vez damos valores específicos para los betas? Esta es una cuestión que se tratará más adelante cuando discutamos los métodos de estimación de parámetros, nombre con el que se conoce al procedimiento estadístico encaminado a elegir los mejores betas para un conjunto de datos observados.
Por lo pronto vamos a proceder de manera más intuitiva, dejando de lado el formalismo matemático, y vamos a ejemplificar un posible modelo ajustado usando valores \(\beta_0=700\) y \(\beta_1 = 1.2\). Hemos elegido estos valores por simple inspección visual así que no esperamos que sean óptimos en ningún sentido, veamos:
# Calculamos las predicciones de "y" usando el modelo
b0 = 700
b1 = 1.2
x= hass_dolar$precio_dolar
y_pred = b0 + b1*x
hass_dolar %>%
ggplot(aes(x= precio_dolar, y= precio_aguacate_kg ))+
geom_point() +
geom_line( aes( x= x, y = y_pred), col = "red")Aunque no estamos seguros de la optimalidad de los betas elegidos, el ajuste basado en este modelo es a simple vista decente, pero no deberíamos depender de esto. Es necesario usar las herramientas que R nos ofrece para estimar los betas. Antes de pasar a su uso ciego, revisemos en qué consideraciones estadísticas se basan estos métodos. Para ello definamos primero los errores cometidos por el modelo, y más aún, la medida de desempeño que usaremos para decidir cuál es la mejor elección de los betas. Usaremos entonces la métrica MSE(Mean Squared Errors) dada por:
\[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}e_i²\]
Por otro lado podemos decir que \(e_i = y_i - \hat{y}\), donde la cantidad \(\hat{y}\) representa el valor de la variable \(Y\) que arroja el modelo, el cual es por regla general diferente al verdadero valor observado de \(Y\), es decir, los errores son las diferencia de lo observado \(y_i\) y lo predicho por el modelo basado en la recta \(\hat{y}\), por lo tanto:
\[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}(y_i - \hat{y})²\]
Y más específicamente:
\[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}(y_i - (\beta_0 + \beta_1*x_i))²\] \[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}(y_i - \beta_0 - \beta_1*x_i)²\]
Con esta expresión ya se hace patente que \(MSE_{\text{NO ingenuo}}\) es función de los pares de valores \((x_i, y_i)\) a los que se debe ajustar la recta, así como de los valores que elijamos para los betas.
Sin embargo, como los pares \((x_i, y_i)\) ya vienen dados por nuestra base de datos o, hablando en términos más generales, por la información recaudada para la construcción del modelo, sólo tendremos opción de variar los betas en aras de minimizar la cantidad \(MSE_{\text{NO ingenuo}}\). Dicho de otra manera los pares de valores \((x_i, y_i)\) se comportan como constantes a lo largo de mi proceso de minimización.
Dado lo anterior, es buena idea definir una función MSE en R que permita calcular, para diferentes elecciones de los betas, la medida de desempeño. Esta medida de desempeño es conocida también en la literatura como función de pérdida, queriendo indicar con ello una pérdida de ajuste del modelo a los datos. Es importante aclarar que no existe una única elección para la función de pérdida de un modelo predictivo, otras alternativas típicas son la mediana de las desviaciones absolutas (MEDA), la media de los errores absolutos (MAE), entre otros. Sin embargo en regresión lineal es usual usar MSE porque su expresión como cuadrado de errores permite usar teoría de inferencia basada en la distribución F cuando los errores se asumen normales. Veremos esto más a detalle cuando discutamos pruebas de hipótesis sobre el modelo de regresión.
1.10 Una relación importante entre los betas del modelo lineal
Retomando nuestro modelo original
\[y = \beta_{0} + \beta_{1} * x + e \text{ con } e \text{~}N(0, \sigma^2)\]
Tomando valor esperado en ambos lados de la ecuación tenemos que:
\[E(y) = E(\beta_{0} + \beta_{1} * x + e)\] Pero como \(\beta_0\) y \(\beta_1\) son constantes, tenemos que:
\[E(y) = E(\beta_0) + E(\beta_1*x) + E(e)\] Así que:
\[E(y) = \beta_0 +\int_{-\infty}^{\infty}\beta_1 *xf(x)dx + E(e)\] \[E(y) = \beta_0 + \beta_1*\int_{-\infty}^{\infty}xf(x)dx + E(e)\] Y como \(\int_{-\infty}^{\infty}xf(x)dx\) es por definición \(E(x)\) y además \(E(e) = 0\) por diseño (de acuerdo con el supuesto usado para la distribución del error), entonces:
\[E(y) = \beta_0 + \beta_1*E(x)\] Dicho de otra manera, la recta del modelo debe pasar por el punto \((E(x), E(y))\) !!
La recta debe pasar por el centroide de la nube de puntos.
Ciertamente no conocemos el valor esperado ni de la variable predictora X ni de la variable respuesta Y, pero buenos estimados basados en información muestral serían \(\bar{x}\) y \(\bar{y}\). Es decir, no podremos saber con precisión los valores \(\beta_0\) y \(\beta_1\), pero nos podemos conformar con buenas estimaciones que cumplan la condición:
\[\bar{y} = \color{red}{\widetilde\beta_0} + \color{red}{\widetilde\beta_1}\bar{x}\]
El coloreado y la tilde ancha para los betas pretenden enfatizar que estos valores son estimados de los correspondientes \(\beta_0\) Y \(\beta_1\) del modelo teórico subyacente. En conclusión: la mejor recta estimada deberá pasar por \((\bar{x}, \bar{y})\), y por lo tanto hay una relación atando a los betas estimados:
\[\begin{equation} \label{eq:ec0} \color{red}{\widetilde\beta_0} = \bar{y} - \color{red}{\widetilde\beta_1}\bar{x} \end{equation} \] Es decir que dando valores a \(\color{red}{\widetilde\beta_1}\) podremos obtener valores para \(\color{red}{\widetilde\beta_0}\) y en últimas dibujar diversas rectas en el plano, con el fin de establecer cuál de ellas minimiza el \(MSE_{\text{NO ingenuo}}\).
1.11 Una animación que ilustra la estimación de betas
A continuación presentamos un código en R que permite comprender cómo funcionan numéricamente los métodos de estimación de betas para un modelo lineal simple, específicamente cuando se usa como método de estimación el criterio de minimización del promedio de los cuadrados de los errores:
# Extraemos las variables de interés
x = hass_dolar$precio_dolar
y = hass_dolar$precio_aguacate_kg
# Calculamos medias muestrales para las variables X y Y
x_barra <- mean(x)
y_barra <- mean(y)
# Definimos posibles valores para b1, y por consiguiente para b0
b1_values <- seq(-1, 3, by = 0.08)
b0_values <- y_barra - b1_values * x_barra
# Fijamos límites en X y Y para cada gráfico
xmin = min(x) - 50
xmax = max(x) + 50
ymin = min(y) - 50
ymax = max(y) + 50
# Construimos una función para calcular el MSE en función de x_i, y_i, b1 y b0
MSE <- function(x, y, b0, b1) {
errors <- y - (b0 + b1*x)
squared_errors <- errors^2
mse <- mean(squared_errors)
return(mse)
}
# Crear un dataframe para almacenar los resultados
list_mse <- vector()
# Calcular el MSE para diferentes valores de B1
for (i in 1:length(b1_values)) {
mse <- MSE(x, y, b0_values[i], b1_values[i])
list_mse[i] <- mse
}
mse_df = data.frame(
B0 = b0_values,
B1 = b1_values,
mse = list_mse,
color_point = 'green'
)
# Capturando el renglón donde ocurre el mínimo mse
pos_min <- which(mse_df$mse==min(mse_df$mse))
# Impriendo las filas cercanas al mínimo
mse_df$color_point[(pos_min - 4):(pos_min + 4)] <- rep('red',9)
ftable(mse_df[(pos_min - 8):(pos_min + 8), ])B0 | B1 | mse | color_point |
|---|---|---|---|
3,475.69602 | 0.28 | 584,101.9 | green |
3,230.39243 | 0.36 | 516,278.7 | green |
2,985.08884 | 0.44 | 457,137.2 | green |
2,739.78525 | 0.52 | 406,677.3 | green |
2,494.48166 | 0.60 | 364,899.1 | red |
2,249.17807 | 0.68 | 331,802.7 | red |
2,003.87448 | 0.76 | 307,387.9 | red |
1,758.57089 | 0.84 | 291,654.8 | red |
1,513.26730 | 0.92 | 284,603.3 | red |
1,267.96371 | 1.00 | 286,233.6 | red |
1,022.66012 | 1.08 | 296,545.5 | red |
777.35653 | 1.16 | 315,539.2 | red |
532.05294 | 1.24 | 343,214.5 | red |
286.74934 | 1.32 | 379,571.5 | green |
41.44575 | 1.40 | 424,610.2 | green |
-203.85784 | 1.48 | 478,330.6 | green |
-449.16143 | 1.56 | 540,732.6 | green |
¿Qué hemos logrado hasta acá?
- Demostrar que el valor de \(MSE_{\text{NO ingenuo}}\) depende de los betas elegidos para el modelo.
- Demostrar cómo el valor de \(\color{red}{\widetilde\beta_0}\) está atado con el de \(\color{red}{\widetilde\beta_1}\) si de verdad queremos que la recta estimada respete la condición de que \(e \text{~}N(0, \sigma^2)\) y por lo tanto de que \(E(e)=0\).
- Generar un dataframe que calcula los valores de \(MSE_{\text{NO ingenuo}}\) para cada par de valores elegidos para \(\color{red}{\widetilde\beta_0}\) y \(\color{red}{\widetilde\beta_1}\), dando valores a \(\color{red}{\widetilde\beta_1}\)
- Estimamos que el \(\color{red}{\widetilde\beta_1}\) óptimo está cerca del valor 0.94, y que el \(\color{red}{\widetilde\beta_0}\) óptimo está cerca del valor 1452.
Ahora graficaremos una curva que muestre esta relación:
if (!file.exists("beta_vs_mse.gif")){
# Crear la animación
animation <- ggplot(mse_df, aes(B1, mse)) +
geom_line(aes( x= B1, y = mse, color = color_point), linetype = "dashed") +
geom_point(aes(color = color_point), size = 3) +
labs(title = "Valores de MSE para cada posible B1, respetando E(e) = 0") +
geom_text(data = mse_df,
aes(label = paste("MSE: ", round(mse,0))), x= 1, y = 2e6, color='black'
) +
geom_text(data = mse_df,
aes(label = paste("MSE mínimo: ", round(min(mse),0))),
x= 1, y = min(mse_df$mse) - 1e5, color='red') +
scale_color_manual(values = c("red" = "red", "green" = "green"),
labels = c("Subóptima", "Óptima"),
name = "Soluciones") +
theme_minimal() +
transition_reveal(B1)
animate(animation, duration = 10, fps = 2, renderer = gifski_renderer())
anim_save("beta_vs_mse.gif")
}Ahora veamos el efecto que tiene la elección de los betas en la recta regresora, e incluyamos pausas en la animación para detenernos en el momento en que encontremos la recta óptima.
if (!file.exists("beta_estimation.gif")){
# Las pausas en la animación son fotogramas con información repetida correspondiente
# al renglón óptimo
fotogramas_pausa <- data.frame(
B0 = rep(b0_values[pos_min], 2*nrow(mse_df)),
B1 = rep(b1_values[pos_min], 2*nrow(mse_df)),
mse = rep(list_mse[pos_min], 2*nrow(mse_df)),
color_point = rep('red',2*nrow(mse_df))
)
mse_df <- rbind(mse_df, fotogramas_pausa)
mse_df <- mse_df[order(mse_df$B1), ]
mse_df$estado <- seq(1,nrow(mse_df))
# Construimos el dataframe de predicciones para cada estado, con el cual poder
# graficar los residuales en cada fotograma
df_predicciones <- data.frame(B0_pred=numeric(0),
B1_pred=numeric(0),
x_from_pred=numeric(0),
y_from_pred=numeric(0),
y_pred_mod=numeric(0),
estado = numeric(0))
for (i in 1:nrow(mse_df)){
df_nuevo <-data.frame(
BO_pred = rep(mse_df[i, 'B0'], length(x)),
B1_pred = rep(mse_df[i, 'B1'], length(x)),
x_from_pred = x,
y_from_pred = y,
y_pred_mod = mse_df[i, 'B0'] + mse_df[i,'B1']*x,
estado = rep(mse_df[i, 'estado'], length(x))
)
df_predicciones <- rbind(df_predicciones, df_nuevo)
}
# Creamos el gráfico base
base_plot <- ggplot(data.frame(x = x, y = y), aes(x = x, y = y)) +
geom_point() +
geom_point(x=mean(x), y=mean(y), color='green', size=2.5) +
geom_abline(aes(slope = 0, intercept = y_barra),
linetype = "dashed", color = "blue") +
labs(title = "Estimación de Betas por Minimización del MSE") +
theme_minimal()
# Creamos la animación
animation <- base_plot +
geom_abline(aes(slope = B1, intercept = B0, color = color_point), data=mse_df) +
geom_text(aes(label = paste("MSE: ", round(mse,0))), x= 2000, y = 6000, data=mse_df,
color = "red")+
geom_segment(data = df_predicciones,
aes(x=x_from_pred, y=y_from_pred, xend=x_from_pred, yend=y_pred_mod),
col = "violet", lty='dashed') +
transition_states(estado, transition_length = 2, state_length = 1) +
scale_color_manual(values = c("red" = "red", "green" = "green"),
labels = c("Subóptima", "Óptima"),
name = "Rectas regresoras") +
enter_fade() +
exit_fade()
# Guardamos la animación en un archivo GIF
animate(animation, duration = 20, fps = 2, renderer = gifski_renderer())
anim_save("beta_estimation.gif", animation)
}1.12 Usando funciones en R para estimar parámetros
Finalmente hemos llegado al punto en que usaremos las funciones de R para estimar parámetros del modelo de regresión lineal. Esto se hará bajo el enfoque de optimización, para ello conviene darle una mirada al siguiente problema:
Estime un par de valores \((\widetilde\beta_0,\widetilde\beta_1)\) tal que minimice la expresión \(MSE = f(\beta_0,\beta_1)\), donde:
\[f(\beta_0,\beta_1) = \frac{1}{n}\sum^{n}_{i=1}(y_i - \beta_0 - \beta_1*x_i)²\] Tomando en cuenta que no hay restricciones aplicadas a los betas.
Para resolver el problema primero efectuemos algunas simplificaciones tomando partida de la siguiente identidad algebraica: \((a+b+c)² = a²+b²+c²+2ab+2ac+2bc\), de tal manera que:
\[f(\beta_0,\beta_1) = \frac{1}{n}\sum^{n}_{i=1}(y_i²+\beta_0²+\beta_1²x_i²-2y_i\beta_0-2y_i\beta_1x_i+2\beta_0\beta_1x_i)\]
Ahora distribuyamos la suma y separemos aparte los términos que dependen de su índice:
\[f(\beta_0,\beta_1) = \frac{1}{n}\left[\sum^{n}_{i=1}y_i²+\beta_0²\sum^{n}_{i=1}1+\beta_1²\sum^{n}_{i=1}x_i²-2\beta_0\sum^{n}_{i=1}y_i-2\beta_1\sum^{n}_{i=1}y_ix_i+2\beta_0\beta_1\sum^{n}_{i=1}x_i\right]\] Simplificando términos y distribuyendo la fracción:
\[f(\beta_0,\beta_1) = \frac{1}{n}\sum^{n}_{i=1}y_i²+\frac{\beta_0²}{n}n+\frac{\beta_1²}{n}\sum^{n}_{i=1}x_i²-\frac{2\beta_0}{n}\sum^{n}_{i=1}y_i-\frac{2\beta_1}{n}\sum^{n}_{i=1}y_ix_i+\frac{2\beta_0\beta_1}{n}\sum^{n}_{i=1}x_i\] Para finalmente obtener:
\[f(\beta_0,\beta_1) = \beta_0² +\bar{x²}\beta_1² -2\bar{y}\beta_0-\left[\frac{2}{n}\sum^{n}_{i=1}x_iy_i\right]\beta_1+2\bar{x}\beta_0\beta_1+\bar{y²} \] En la expresión anterior, las cantidades \(\bar{x²}\), \(-2\bar{y}\), \(-\frac{2}{n}\sum^{n}_{i=1}x_iy_i\), \(2\bar{x}\), \(\bar{y²}\) son todas constantes por depender sólo de los datos que alimentan al modelo, en ese sentido son coeficientes numéricos para los respectivos betas de la función.
LLegados a este punto, podemos sacar las siguientes conclusiones:
La función \(f\) sólo depende de los betas, puesto que los demás valores son constantes, es decir los valores \((x_i, y_i)\) son números consignados en una base de datos, y \(n\) es la cantidad de pares de números consignados. Todos los coeficientes dependerán de operaciones sobre estos números.
De acuerdo con lo anterior la función \(f\) es una función multivariada, pues depende de dos variables: \(\beta_0\) y \(\beta_1\).
Existen métodos muy sencillos para resolver un problema de optimización no restringido cuando la función a optimizar es convexa, entonces sólo nos resta demostrar que la función anterior lo es.
Definición de convexidad
Para demostrar que \(f(\beta_0, \beta_1)\) es convexa, debemos verificar que su matriz hessiana sea semidefinida positiva. La matriz hessiana es una matriz de segundas derivadas parciales de la función.
Primero, calculemos las derivadas parciales de \(f\) con respecto a \(\beta_0\) y \(\beta_1\):
\[\begin{equation} \label{eq:ec1} \frac{\partial f}{\partial \beta_0} = 2\beta_0-2\bar{y}+2\bar{x}\beta_1 \end{equation} \] \[\begin{equation} \label{eq:ec2} \frac{\partial f}{\partial \beta_1} = 2\bar{x²}\beta_1-\frac{2}{n}\sum^{n}_{i=1}x_iy_i+2\bar{x}\beta_0 \end{equation} \] Luego, calculemos las segundas derivadas parciales de \(f\) con respecto a \(\beta_0\) y \(\beta_1\):
\[\frac{\partial^2 f}{\partial \beta_0^2} = 2 > 0\] \[\frac{\partial^2 f}{\partial \beta_1^2} = 2\bar{x^2} > 0\] \[\frac{\partial^2 f}{\partial \beta_0 \partial \beta_1} = 2\bar{x}\]
La matriz hessiana de \(f\) es:
\[ H(f) = \begin{bmatrix} \frac{\partial^2 f}{\partial \beta_0^2} & \frac{\partial^2 f}{\partial \beta_0 \partial \beta_1} \\ \frac{\partial^2 f}{\partial \beta_0 \partial \beta_1} & \frac{\partial^2 f}{\partial \beta_1^2} \end{bmatrix} \]
Por lo tanto:
\[ H(f) = \begin{bmatrix} 2 & 2\bar{x}\\ 2\bar{x} & 2\bar{x^2} \end{bmatrix} \]
Para que \(f(\beta_0, \beta_1)\) sea convexa, la matriz hessiana debe ser semidefinida positiva, lo que significa que debe cumplir la propiedad: \[v^ \top Hv>=0 \text{ para }\forall v \in \mathbb{R}² \] Es decir, debemos verificar si:
\[\begin{bmatrix} a & b \end{bmatrix}H\begin{bmatrix} a \\ b \\ \end{bmatrix}>=0 \text{ para } \forall (a,b) \]
Veamos si esto se cumple:
\[ \begin{bmatrix} a & b \end{bmatrix} \begin{bmatrix} 2 & 2\bar{x} \\ 2\bar{x} & 2\bar{x^2} \end{bmatrix} \begin{bmatrix} a \\ b \\ \end{bmatrix}\\= \begin{bmatrix} 2a+2b\bar{x} & 2a\bar{x} + 2b\bar{x²} \end{bmatrix} \begin{bmatrix} a \\ b \\ \end{bmatrix}\\=2a²+2ab\bar{x}+2ab\bar{x}+2b²\bar{x²} \\=2a²+4ab\bar{x}+2b²\bar{x²}\\ =2(a²+2ab\bar{x}+b²(\bar{x})²)-2b^2(\bar{x})²+2b²\bar{x²}\\ =2(a+b\bar{x})²+2b²(\bar{x²}-(\bar{x})²) \]
Acá es interesante notar que:
\[\bar{x²}-(\bar{x})²=\frac{1}{n}\left[\sum^{n}_{i=1}x_i²\right]-\bar{x}\bar{x} =\frac{1}{n}\left[\sum^{n}_{i=1}x_i²-\bar{x}(n\bar{x})\right]\\ =\frac{1}{n}\left[\sum^{n}_{i=1}x_i²-\bar{x}\sum^{n}_{i=1}x_i\right] =\frac{1}{n}\left[\sum^{n}_{i=1}(x_i²-\bar{x}x_i)\right] \\ =\frac{1}{n}\left[\sum^{n}_{i=1}((x_i²-2\bar{x}x_i)+\bar{x}x_i)\right] =\frac{1}{n}\left[\sum^{n}_{i=1}\left((x_i²-2\bar{x}x_i+\bar{x}²)+(\bar{x}x_i-\bar{x}²)\right)\right]\\ =\frac{1}{n}\left[\sum^{n}_{i=1}\left((x_i-\bar{x})²+\bar{x}(x_i-\bar{x})\right)\right] =\frac{1}{n}\sum^{n}_{i=1}(x_i-\bar{x})²+\bar{x}\left[\sum^{n}_{i=1}\frac{x_i}{n}-\sum^{n}_{i=1}\frac{\bar{x}}{n}\right]\\ =\sigma²_{x}+\bar{x}\left(\bar{x}-\bar{x}\sum^{n}_{i=1}\frac{1}{n}\right)\\ =\sigma²_{x}+\bar{x}(\bar{x}-\bar{x}*1)=\sigma²_{x}\]
Obteniendo finalmente
\[\begin{equation} \label{eq:sigma_x} \bar{x²}-(\bar{x})² = \sigma²_{x} \end{equation} \]
Luego:
\[v^ \top Hv = 2(a+b\bar{x})²+2b²\sigma²_{x} \text{ con } \text{ }v=(a,b)\] Y como ambos sumandos son no negativos se tiene que \[v^ \top Hv >= 0 \text{ para }\forall v\in \mathbb{R}²\] Es decir, la matriz hessiana es semidefinida positiva y por lo tanto \(f(\beta_0,\beta_1)\) es convexa en el espacio de parámetros \(\beta_0\) y \(\beta_1\). Esto implica que el problema de regresión lineal, que busca minimizar esta función de errores, tiene un mínimo global y puede ser resuelto de manera eficiente mediante métodos de optimización convexa.
Una condición necesaria y suficiente para que \((\widetilde\beta_0,\widetilde\beta_1)\) sea un mínimo global es que \(\nabla{} f(\widetilde\beta_0,\widetilde\beta_1)=0\)
Por lo tanto retomando \(\ref{eq:ec1}\) y \(\ref{eq:ec2}\) el problema de minimización queda resuelto hallando la solución al sistema de ecuaciones:
\[\begin{equation} \label{eq:ec3} 2\widetilde\beta_0-2\bar{y}+2\bar{x}\widetilde\beta_1=0 \end{equation} \] \[\begin{equation} \label{eq:ec4} 2\bar{x²}\widetilde\beta_1-\frac{2}{n}\sum^{n}_{i=1}x_iy_i+2\bar{x}\widetilde\beta_0=0 \end{equation} \]
De modo que despejando \(\widetilde\beta_0\) en \(\ref{eq:ec3}\) tenemos lo siguiente:
\[\begin{equation} \label{eq:ec5} \widetilde\beta_0 = \bar{y} - \widetilde\beta_1\bar{x} \end{equation} \] Esto en concordancia con \(\ref{eq:ec0}\), y además haciendo reemplazos en \(\ref{eq:ec4}\) obtenemos:
\[ 2\bar{x²}\widetilde\beta_1-\frac{2}{n}\sum^{n}_{i=1}x_iy_i+2\bar{x}(\bar{y}-\widetilde\beta_1\bar{x})=0\\ 2\bar{x²}\widetilde\beta_1-2\bar{x}²\widetilde\beta_1 = \frac{2}{n}\sum^{n}_{i=1}x_iy_i-2\bar{x}\bar{y}\\ (\bar{x²} - \bar{x}²)\widetilde\beta_1 = \frac{1}{n}\sum^{n}_{i=1}x_iy_i-\bar{x}\bar{y} \] Ahora, usando el hecho de que el lado derecho de la ecuación anterior es por definición \(cov(x,y)\), es decir, la covarianza entre la variable regresora \(x\) y la variable respuesta \(y\), así como la expresión dada en \(\ref{eq:sigma_x}\) tenemos:
\[\begin{equation} \label{eq:ec6} \widetilde\beta_1 = \frac{cov(x,y)}{\sigma_x²} \end{equation} \] Las ecuaciones \(\ref{eq:ec5}\) y \(\ref{eq:ec6}\) son muy informativas en su estructura, la primera indica que la recta regresora debe pasar por el centroide de la nube de puntos: \((\bar{x}, \bar{y})\) y la segunda indica que la pendiente de la recta es proporcional al grado de relación lineal que tenga la variable \(x\) con la variable \(y\) medida a través de la covarianza.
Estos sencillos cálculos son realizados por R a través de la función “lm”, así:
## [1] 0.9449775
## [1] 1436.679
## [1] 288485.9
# Confirmamos usando la función lm, a continuación la documentación
# lm(formula, data, subset, weights, na.action,
# method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
# singular.ok = TRUE, contrasts = NULL, offset, ...)
mod1 <- hass_dolar %>% lm(precio_aguacate_kg ~ precio_dolar,.)
print(mod1$coefficients)## (Intercept) precio_dolar
## 1436.6790256 0.9449775
## [1] 537.1089
parametros <- data.frame(metodo = c("manual","lm"),
beta_0 = c(b0,mod1$coefficients[1]),
beta_1 = c(b1,mod1$coefficients[2]),
sigma_2 = c(sigma_2, (resumen$sigma)^2))
parametros %>%
ftable()metodo | beta_0 | beta_1 | sigma_2 |
|---|---|---|---|
manual | 1,436.679 | 0.9449775 | 288,485.9 |
lm | 1,436.679 | 0.9449775 | 288,485.9 |
#print(mod1$)
# Generando el modelo con R Base sería así:
# mod1 <- lm(promedio_aguacate_kg ~ promedio_dolar, hass_dolar)Podemos notar entonces que ambos métodos son consistentes entre si, algo que no es de extrañar pues internamente R implementa estos métodos en sus librerías.
1.13 Representación funcional del modelo y comprobación de supuesto de normalidad.
Hemos llegado al punto en que ya tenemos todos los elementos necesarios para definir formalmente el modelo estimado sobre nuestros datos. Podemos ahora concluir que el mejor modelo para describir la relación entre la variable precio del aguacate (\(y\)) y precio del dólar(\(x\)), pensado desde la perspectiva de minimización de cuadrados del error y estando restringidos a la familia de modelos de regresión lineal con componente de error normal es:
\(y = 1436.68 + 0.94x + e \text{ con e~N(0,288485.9)}\)
Pero, aunque nos hemos esforzado por encontrar las mejores estimaciones de los betas para minimizar la función de pérdida, esto no garantiza que nuestro modelo induzca errores normales, por tal motivo conviene recordar lo siguiente:
💭 Cuando planteamos el modelo, tuvimos que suponer algunas cosas, entre ellas, que la distribución de los errores sería normal. Por lo tanto, ahora que tenemos construído nuestro modelo, es necesario que validemos la suposición hecha.
🤔💭 ¿Cómo hacerlo ❔
Consideremos la siguiente tabla comparativa que informa sobre las ventajas y desventajas de las tres pruebas de normalidad (Kolmogorov-Smirnov, Lilliefors y Shapiro-Wilk) más usadas en estadística
tabla_normalidad <- data.frame(
"Prueba" = c("Kolmogorov-Smirnov (KS)", "Prueba de Lilliefors", "Prueba de Shapiro-Wilk (SW)"),
"Ventajas" = c(
"- Puede utilizarse para cualquier distribución teórica.\n- Es una prueba no paramétrica versátil.\n - Adecuada para muestras grandes.",
"- Es específica para evaluar la normalidad.\n- Puede ser más adecuada para muestras pequeñas.\n - Utiliza estadísticas de prueba que se ajustan a la distribución estándar.",
"- Es especialmente potente para muestras pequeñas.\n- Diseñada específicamente para evaluar la normalidad.\n- Sensible a desviaciones de la normalidad en cualquier parte de la distribución."
), "Desventajas" = c(
"- Puede ser menos poderosa en muestras pequeñas.\n - No es específica para la distribución normal.",
"- Restringida a la comprobación de normalidad. \n - Menos potente en muestras grandes.\n- No es tan versátil como KS para otras distribuciones.",
"- Más compleja de entender y aplicar.\n- No proporciona estimaciones de parámetros.\n- Puede ser menos adecuada para muestras muy grandes.\n- No es tan versátil como KS para otras distribuciones."))
# Imprimir la tabla con formato
tabla_normalidad %>%
ftable() %>%
align(i = NULL, j = c(2:3),
align = "left", part = "all")Prueba | Ventajas | Desventajas |
|---|---|---|
Kolmogorov-Smirnov (KS) | - Puede utilizarse para cualquier distribución teórica. | - Puede ser menos poderosa en muestras pequeñas. |
Prueba de Lilliefors | - Es específica para evaluar la normalidad. | - Restringida a la comprobación de normalidad. |
Prueba de Shapiro-Wilk (SW) | - Es especialmente potente para muestras pequeñas. | - Más compleja de entender y aplicar. |
Recordemos que la elección de la prueba depende de varios factores, incluyendo el tamaño de la muestra, el conocimiento previo sobre la distribución de los datos y los objetivos del análisis. Ninguna prueba es perfecta y es importante considerar el contexto y las características específicas de tus datos al seleccionar una prueba de normalidad.
🤔💭 Y si el texto anterior dice: “el conocimiento previo sobre la distribución de los datos” 🤔💭 ¿Será necesario representar la distribución en un gráfico? ¿Cómo funciona la prueba Lilliefors? A continuación describiremos los cálculos implicados:
# Paso 1: Ordenamos los datos de menor a mayor
n <- length(mod1$residuals)
df_residuals <- data.frame(residuals = mod1$residuals)
df_residuals <- df_residuals %>%
arrange(residuals)
# Paso 2: Estimamos los parámetros para la distribución normal teórica a partir de los datos
mean_residuals <- mean(df_residuals$residuals)
sd_residuals <- sd(df_residuals$residuals)
# Paso 3: Calculamos la función acumulada para cada valor de residual
cdf_teorica <- pnorm(df_residuals$residuals,
mean = mean_residuals,
sd = sd_residuals)
# Paso 4: Calculamos la máxima diferencia absoluta entre la densidad empírica y la teórica
# Método 1: Construir la función empírica usando fórmulas de R (ecdf) y luego tomar diferencias
# absolutas entre valor empírico y teórico (enfoque KS)
cdf_empirica <- ecdf(df_residuals$residuals)
diferencias_absolutas <- abs(cdf_empirica(df_residuals$residuals) - cdf_teorica)
pos_estadistico_m1 <- which.max(diferencias_absolutas)
estadistico_m1 <- diferencias_absolutas[pos_estadistico_m1]
print(paste0("Este es el valor del estadístico: ", estadistico_m1," y esta es la posición en ",
"la que ocurre la máxima diferencia ", pos_estadistico_m1))## [1] "Este es el valor del estadístico: 0.0566890639519794 y esta es la posición en la que ocurre la máxima diferencia 60"
# Método 2: Construir manualmente los cálculos que permiten hallar la máxima diferencia
# entre la empírica y la teórica usando la definición de la empírica como
# función a trozos 📈
e_plus <- seq(1:n)/n
e_minus <- seq(0:(n-1))/n
D_plus <- e_plus - cdf_teorica
D_minus <- cdf_teorica - e_minus
pos_max_d_plus <- which.max(D_plus)
pos_max_d_minus <- which.max(D_minus)
estadistico_m2 <- max(D_plus[pos_max_d_plus], D_minus[pos_max_d_minus])
pos_estadistico_m2 <- function() {
if(max(D_plus)>=max(D_minus)) {
return(pos_max_d_plus)
}else{
return(pos_max_d_minus)
}
}
print(paste0("Este es el valor del estadístico: ", estadistico_m2," y esta es la posición en ",
"la que ocurre la máxima diferencia ", pos_estadistico_m2()))## [1] "Este es el valor del estadístico: 0.0566890639519794 y esta es la posición en la que ocurre la máxima diferencia 60"
# Comprobamos los resultados invocando la función de R que realiza estos cálculos automáticos
Lilliefors_test <- lillie.test(df_residuals$residuals)
print(Lilliefors_test)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: df_residuals$residuals
## D = 0.056689, p-value = 0.3637
estadistico_l <- Lilliefors_test$statistic
print(paste0("Este es el estadístico de Lilliefors usando paquetería de R: ", estadistico_l))## [1] "Este es el estadístico de Lilliefors usando paquetería de R: 0.0566890639519794"
Hasta acá hemos verificado que los cálculos realizados manualmente coinciden con el código implementado por la librería nortest de R, y más específicamente en el código fuente de la función lillie.test. Ahora ilustremos algo de estos pasos usando una gráfica:
x_seg <- df_residuals$residuals[pos_estadistico_m1]
y_sup_seg <- cdf_teorica[pos_estadistico_m1]
y_inf_seg <- cdf_empirica(x_seg)
cat("x_seg: ", format(x_seg, digits=5),
" y_inf_seg: ", format(y_inf_seg, digits=5),
" y_sup_seg: ", format(y_sup_seg, digits=5),
sep='')## x_seg: -147.97 y_inf_seg: 0.44776 y_sup_seg: 0.39107
datos_para_grafico <- data.frame(
residuales = df_residuals$residuals
)
ggplot(datos_para_grafico, aes(x = residuales)) +
# geom_density(color = "blue", fill = "lightblue", alpha = 0.5) +
stat_function(fun = pnorm, args = list(mean = mean_residuals,
sd = sd_residuals),
color = "red", size = 1) +
stat_ecdf(color = "green", size = 1) +
geom_segment(aes(x = x_seg, y = y_inf_seg, xend = x_seg,
yend = y_sup_seg),
color = "purple", size = 1, linetype = "solid") +
geom_point(aes(x = x_seg, y = y_inf_seg),
color = "purple",size = 3) +
geom_point(aes(x = x_seg, y = y_sup_seg),
color = "purple", size = 3) +
annotate("text", x = x_seg , y = (y_inf_seg + y_sup_seg)/2,
label = paste0("Máxima Separación=",
format(estadistico_m1,
digits=5)," 🤔💭", sep=''),
color = "purple", size = 3, hjust = -0.05, vjust = 0) +
coord_cartesian(xlim = c(min(datos_para_grafico$residuales),
max(datos_para_grafico$residuales))) +
xlab("Residuales del modelo") +
ylab("Probabilidad acumulada") +
labs(title = "Esquema de funcionamiento de la prueba de Lilliefors") +
theme_minimal()Nos queda la siguiente reflexión sobre la gráfica anterior, ¿Por qué elegir la máxima diferencia absoluta entre las dos curvas en lugar de un promedio de las diferencias a lo largo de todo el eje? ¿Por qué no usar una media ponderada de los residuales al cuadrado?
Comparemos los resultados contra la prueba de Shapiro WIlk y Kolmogorov
##
## Shapiro-Wilk normality test
##
## data: mod1$residuals
## W = 0.98361, p-value = 0.1078
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: mod1$residuals
## D = 0.50746, p-value < 0.00000000000000022
## alternative hypothesis: two-sided
Hemos notado diferencias con la prueba de Kolmogorov, ¿Cuál es la razón?, podrás encontrar en el siguiente artículo en el que Lilliefors introduce por primera vez su prueba:
Lilliefors, H. (June 1967), “On the Kolmogorov–Smirnov test for normality with mean and variance unknown”, Journal of the American Statistical Association, Vol. 62. pp. 399–402.
Una aclaración respecto a que esta prueba es menos potente ❗ en comparación con la prueba KS estándar, esto significa que es menos sensible para detectar desviaciones de la normalidad en muestras grandes, especialmente si los datos no siguen una distribución estrictamente normal. La prueba KS estándar es más robusta en este sentido, pero por lo mismo más conservador. En el contexto de modelos de regresión, no requerimos un cumplimiento estricto de la normalidad sino una semejanza razonable para habilitar el uso del modelo. En ese sentido es preferible usar Lilliefors.😀😎
🤔💭 Otras posibles preguntas son:
- ¿Cuáles son los riesgos que implica asumir a priori un comportamiento específico para el error?
- ¿Si el modelo no cumple con lo supuesto para el error, qué alternativas tengo?
- Al plantear un modelo alternativo que captura mejor el comportamiento del error usando una distribución distinta a la normal ¿Puedo seguir usando los métodos de inferencia tradicionales o qué modificaciones debo introducir?
- ¿Por qué es importante para un magister en estadística conocer sobre métodos de simulación Montecarlo?
2. EL SUPUESTO DE INDEPENDENCIA DEL ERROR
Un supuesto que hemos perdido de vista por un momento en los apartados anteriores es el supuesto de independencia del error. Su consideración es importante porque representa un punto de quiebre en el uso de modelos de regresión vs modelos alternativos como los de series de tiempo. Pero ¿Qué significa exactamente la independencia?. Consideremos los siguientes escenarios:
Un modelo que presenta sistemáticamente más dificultades para predecir de forma exacta valores altos de la variable respuesta pero es bueno prediciendo valores bajos.
Un modelo que se construyó sobre datos de un experimento de durezas de metales, donde el mecanismo que sensa la dureza perdió calibración o exactitud después de tomar la medida de las 100 primeras muestras.
Un modelo que intento predecir un precio de activo financiero en función de un activo de referencia, pero conforme pasa el tiempo el modelo se va desajustando y haciendo menos preciso.
Un modelo que desconce la presencia de una variable oculta que controla la variación acoplada de mi variable respuesta vs mi variable predictora y las hace parecer más correlacionadas de lo que realmente son (regresiones espurias)
Para cada uno de los casos anteriores el investigador tiene que decidir el tipo de recurso gráfico o estadístco que le permita inspeccionar si una situación de esa naturaleza está ocurriendo. Las alternativas son:
Gráficos de residuales vs. valores ajustados: Graficar los residuales (diferencia entre los valores observados y los valores predichos) contra los valores ajustados (las predicciones del modelo) es una forma común de detectar patrones de no independencia.
🤔💭 ¿Qué comportamiento es deseable para este gráfico?
Si no hay patrones evidentes en el gráfico y los residuales se distribuyen aleatoriamente alrededor de cero, es un indicio de independencia.
Gráfico de residuales con índice temporal: Si tus datos están ordenados en el tiempo (series temporales), un gráfico de residuales en el tiempo puede revelar patrones de autocorrelación. Deben parecer ruido blanco, es decir, sin patrones visibles. También es conveniente usar gráficos de autocorrelación de residuales (ACF) y gráficos de autocorrelación parcial de residuales (PACF), ya que estos informan sobre una posible correlación entre el error actualmente observado y errores observados en otros momentos del pasado.
Prueba de Durbin-Watson: Esta prueba estadística evalúa si existe autocorrelación de primer orden en los residuales, mejora el análisis gráfico por ser una prueba formal. Un valor de Durbin-Watson cercano a 2 indica independencia, mientras que valores significativamente diferentes de 2 pueden indicar autocorrelación. De nuevo es necesario elegir una variable de ordenamiento para los errores encaminada a detectar el tipo de patrón de interés.
🤔💭 ¿Cómo hacer esto?
- Ordene los errores según magnitud del valor predicho (ajustado)
- Ordene los errores según el tiempo.
- Ordene los errores según el orden de experimentación.
- No ordene por la posición del registro en el data frame a menos qué entienda qué significa ese orden y qué aporta al análisis.
- Ordene por una variable de interés en la que quiera inspeccionar si los residuales se ven afectados por la magnitud de tal variable.
👦👴👧👨👨👵
Ejemplo: Puede considerar revisar si el consumo calórico puede explicar residuales de un modelo simplificado en el que la estatura pretenda explicar el peso. 🔎🕵️ Aunque este ejemplo parece obvio ilustra una idea muy importante, porque nos da un criterio para ingresar una nueva variable a un modelo de regresión: Ingrésela cuando dicha variable se relaciona bien con los residuales generados por la primera versión de mi modelo!!!. Es decir, la variable debe ingresar al modelo cuando pueda explicar la variabilidad restante no explicada por mi modelo anterior.
Prueba de Ljung-Box: Esta prueba se utiliza para evaluar la autocorrelación en los residuales en varios retardos. Si los residuales son independientes, las autocorrelaciones en los retardos deberían ser cercanas a cero. Por varios retardos entendemos inspeccionar relación no sólo con el valor del residual inmediatamente anterior, sino contra cualquier otro ubicado en un momento más anterior. 🔎 Es una mejora a la prueba de Durbin-Watson en la medida en que no limita la autocorrelación al retardo anterior.
Estas son algunas de las formas comunes de comprobar la independencia de los errores en un modelo de regresión lineal. Si encuentras evidencia de autocorrelación o patrones en los residuales, puede ser necesario considerar modelos más avanzados, como 📉 🚪modelos autorregresivos, o de media móvil, o en forma más general modelos de series temporales, para capturar adecuadamente la estructura de dependencia en los datos.
En conclusión:
En regresión lineal asumimos la no correlación de los errores. En nuestro caso el precio del aguacate y el precio del dólar son variables evolucionando en el tiempo, así que quisiéramos que tras el transcurrir del tiempo nuestro modelo no se deteriore. Una forma en que esto puede ocurrir es que alguna autocorrelación temporal provoque residuales que van variando en sus propiedades conforme el tiempo pasa. Por tal motivo necesitamos técnicas como las anteriores para probar independencia.
Por ejemplo puede ocurrir que el error se vaya volviendo más grande en promedio con el paso del tiempo o vaya creciendo en dispersión, o cualquier otra anomalía que impida suponer que tenemos un modelo estable para todo momento del tiempo. Los gráficos de dispersión típicos ilustran la asociación entre la magnitud de la variable predictora y la variable respuesta, pero para probar independencia es preferible que los errores se dibujen indexados por la fecha en la que tal error ocurrió, y así poder apreciar posibles patrones temporales. Haremos esas inspecciones a continuación 🔎 :
# Primero inspeccionamos las variables individuales para observar su
# comportamiento temporal en un mismo gráfico
# library(lubridate)
x <- hass_dolar$precio_dolar
y <- hass_dolar$precio_aguacate_kg
hass_dolar$fecha <-as.Date(paste0(hass_dolar$anio,"-",hass_dolar$mes,"-01"))
hass_dolar$predicciones <- mod1$fitted.values
hass_dolar$residuales <-mod1$residuals
# Gráfico de evolución de precios y predicciones
ggplot(hass_dolar, aes(x = fecha)) +
geom_line(aes(y = precio_dolar, color = "Precio del dólar"), size = 1) +
geom_line(aes(y = precio_aguacate_kg, color = "Precio del aguacate"), size = 1) +
geom_line(aes(y = predicciones, color = "Predicciones precio aguacate"), size = 1) +
scale_color_manual(values = c("Precio del dólar" = "blue",
"Precio del aguacate" = "green",
"Predicciones precio aguacate" = "red")) +
labs(title = "Evolución temporal del precios y predicciones",
x = "Fecha",
y = "Valor") +
theme_minimal() +
scale_x_date(date_breaks = "1 year", date_labels = "%b %Y") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))De entrada, podemos observar en el gráfico anterior que hay zonas donde el precio del aguacate fue más volátil que el precio del dólar, razón por la cual el precio del dólar no pudo seguir eficientemente el patrón del precio del aguacate, veremos esto cómo afectó los residuales:
# Gráfico de residuales vs valores ajustados
hass_dolar %>%
ggplot(aes(x= predicciones, y=residuales)) +
geom_point()+
theme_light()+
ylab("Residuales") +
xlab("Valor predicho")# Gráfico de residuales con índice temporal
hass_dolar %>%
ggplot(aes(x= fecha, y=residuales)) +
geom_point()+
geom_line(aes(y = residuales, color = "residuales"), size = 1)+
theme_light()+
xlab("Fecha")+
ylab("Residuales")# Prueba de Durbin-Watson tomando el mismo orden en el que venían los
# datos en el dataframe puesto que este orden se corresponde con el
# orden temporal.
resultado_dw <- dwtest(mod1)
print(resultado_dw)##
## Durbin-Watson test
##
## data: mod1
## DW = 0.47441, p-value < 0.00000000000000022
## alternative hypothesis: true autocorrelation is greater than 0
# Prueba de Ljung-Box
# Como esperamos cierta estacionalidad intra año en los datos podría tener sentido esperarla ttambién para los residuales. Por lo tanto eligiremos lag = 1,2,3,...,12
p_values <- vector()
for (i in 1:12) {
resultado_ljung_box <- Box.test(mod1$residuals, lag = i,
type = "Ljung-Box")
p_values[i] <-resultado_ljung_box$p.value
}
resultados_lb <- data.frame(
lag = seq(1,12),
p_value = p_values
)
print(resultado_ljung_box)##
## Box-Ljung test
##
## data: mod1$residuals
## X-squared = 124.84, df = 12, p-value < 0.00000000000000022
## lag p_value
## 1 1 0
## 2 2 0
## 3 3 0
## 4 4 0
## 5 5 0
## 6 6 0
## 7 7 0
## 8 8 0
## 9 9 0
## 10 10 0
## 11 11 0
## 12 12 0
Conclusiones:
Los errores no presentan un patrón claro al momento de compararse contra los valores ajustados, es decir que el modelo tiene una calidad similar prediciendo valores pequeños, intermedios y grandes en el precio del aguacate.
El gráfico de residuales indexados por el tiempo parece tener ciertas oscilaciones predecibles pero es difícil de establecer por simple inspección visual.
La prueba de Durbin Watson nos ayuda a revisar si hay autocorrelación del error de hoy con el error de hace un mes (recuerde que la serie es mensual). De hecho la prueba arroja un estadístico de prueba \(DW < 2\) indicando con ello autocorrelación positiva de los errores.
En consistencia con lo anterior la prueba de Ljung-Box arroja valores p sistemáticamente menores que 0.05 para todos los rezagos considerados: 1 mes, 2 meses, 3 meses, …, 1 año. Y decimos que en consistencia con la prueba de Durbin Watson porque la prueba de Ljung box será indicativa de autocorrelación en máximo n rezagos si ya lo es para máximo n=1 rezagos como lo indica la prueba DW.
La no independencia del error invalida completamente mi modelo, porque en estos escenarios la varianza del error del modelo se encuentra subestimada, pues ante la presencia de autocorrelación de los residuales no podemos afirmar que \(\sigma_e² = var(error)\) ya que esto asume implícitamente independencia. En procesos autocorrelacionados la varianza es realmente mayor que esto. Usar un error subestimado puede llevarme a conclusiones falsamente significativas, como por ejemplo puede declarar la variable predictora como significativa, al estimar inadecuadamente el error para \(\beta_1\) ya que \(\widetilde{V}(\beta_1)\) depende de \(\widetilde{\sigma_e²}\)
Dejaremos el estudio de las cuestiones relativas al tratamiento de errores autocorrelacionados para el capítulo de series de tiempo.
2. EL SUPUESTO DE HOMOCEDASTICIDAD
De manera informal este supuesto busca verificar que la varianza de los errores es estable, bien sea respecto al tiempo, respecto a una variable predictora, o respecto a la magnitud de los valores predichos. Las inspecciones de estas cuestiones usualmente se basan en el estudio de los patrones de los gráficos de residuales, en específico estudiamos:
Si un gráfico de residuales estandarizados (divididos por su desviación estándar) vs valor predicho arroja dispersiones diferentes conforme el valor predicho crece. En un gráfico como éste si observamos que la dispersión de los puntos cambia a lo largo de los valores ajustados (por ejemplo, los puntos se abren o se estrechan a medida que aumentan los valores ajustados), eso podría indicar la presencia de heterocedasticidad.
Si el modelo posee una variable categórica predictiva (no es nuestro caso), respecto a la cual podamos validar si la dispersión del error es mayor en los grupos definidos por esa variable predictiva. Ampliaremos estas cuestiones en el apartado de modelo de regresión lineal múltiple.
Si un gráfico de residuos vs variables independientes revela que la dispersión de los residuos cambia con respecto a cada variable independiente. Por ejemplo, podemos patrones de abanico o embudo que puedan indicar problemas de heterocedasticidad.
Si alguna pruebas estadísticas formales como la prueba de Breusch-Pagan o la prueba de White proporcionan evidencia de heterocedasticidad al comprobarse que los cuadrados de los residuales se relación sistemáticamente con el conjunto de variables predictoras, esto es, si es posible construir un modelo no ingenuo que haga depender el cuadrado del error de los valores de las variables predictoras. Si tal cosa existe significa que los errores son heterocedasticos en la medida en que su varianza depende o cambia según los valores de las variables predictoras. Apliquemos estos análisis a continuación:
# Residuales estandarizados vs valores predichos
hass_dolar %>%
ggplot(aes(x= predicciones, y=residuales/sd(residuales))) +
geom_point()+
theme_light()+
ylab("Residuales estandarizados") +
xlab("Valor predicho")# Residuales estandarizados vs precio del dólar
hass_dolar %>%
ggplot(aes(x= precio_dolar, y=residuales/sd(residuales))) +
geom_point()+
theme_light()+
ylab("Residuales estandarizados") +
ylab("Precio del dólar")Como se puede apreciar, el gráfico de errores estandarizados vs valores predichos será muy parecido al gráfico de errores vs variable predictora en un modelo de regresión lineal simple, puesto que en este contexto las predicciones dependen de la variable de pronóstico en forma tal que la primera no es más que un rescalado y desplazamiento de la segunda. Tendría más sentido en un modelo de regresión lineal múltiple generar ambos gráficos, acá uno sólo de los anteriores bastaría.
Por otro lado, ambos gráficos revelan que no existe un patrón específico que haga sospechar de la presencia de heterocedasticidad, hagamos un inspección más cuidadosa usando pruebas formales:
# Pruebas estadísticas formales
resultado_breusch_pagan <- bptest(mod1)
print(resultado_breusch_pagan)##
## studentized Breusch-Pagan test
##
## data: mod1
## BP = 1.741, df = 1, p-value = 0.187
La prueba de Breusch-Pagan propuesta en (1979) consiste en ajustar un modelo de regresión lineal con variable respuesta dada por residuales del modelo original al cuadrado \(e_i²\) y como covariables las variables del modelo original. Por ejemplo, si se tienen \(k=2\) covariables para explicar a \(Y\) entonces el modelo de regresión para estudiar la homocedasticidad es:
\[e_i² = \delta_0 + \delta_1*x_1 + \delta_2*x_2 + u\]
Si se concluye que \(\delta_1 = \delta_2 = 0\), significa que los residuales no son función de las covariables del modelo. El estadístico en esta prueba está dado por \(n*R²\) y bajo la hipótesis nula verdadera, el estadístico tiene distribución \(\chi²_k\)
Lo que el resultado nos está indicando es que \(\chi²_1 = 1.741\), y \(df = 1\), con \(\text{p value} = 0.187\), y debido a que p>0.05, no hay evidencia estadística suficiente para el rechazo de la hipótesis de homocedasticidad. Todo esto en concordancia con lo revelado por el análisis gráfico.
Sin embargo, el test de Breusch-Pagan sólo detecta formas lineales de heterocedasticidad, esto significa que sólo relaciona los cuadrados del error con términos que dependen linealmente de las covariables. Para resolver esta limitación, el test de White propuesto por White (1980), permite contrastar no linealidades utilizando los cuadrados y los productos cruzados de todos los regresores. Sea que se use Breusch-Pagan o White, note que sólo se están explorando dos formas de heterocedasticidad en un universo infinitamente amplio de formas en que ésta se puede presentar, esto es, sólo se explora heterocedasticidad con forma parabólica y con forma lineal. Se podrían explorar otras pero conforme la complejidad del modelo crece se pierde su utilidad en la práctica.
En el caso del test de White, si \(k=2\) el modelo que relaciona \(e_i²\) con las covariables queda así:
\[e_i² = \delta_0 + \delta_1x_1 + \delta_2x_2 + \delta_3x_1x_2 + \delta_4 x_1² + \delta_5x_2² + u\] Podemos pensar el modelo de White como una extensión del modelo de Breusch-Pagan, por eso no es de extrañar que en la siguiente celda usemos la misma función en R, pero agregando un término al cuadrado para la única covariable implicada en nuestro ejemplo. No habrán términos cruzados porque sólo tenemos una variable predictora.
resultado_white <- bptest(mod1,
varformula = ~ (precio_dolar +
I(precio_dolar^2)),
data=hass_dolar)
print(resultado_white)##
## studentized Breusch-Pagan test
##
## data: mod1
## BP = 7.8925, df = 2, p-value = 0.01933
Para ilustrar como trabajan estas pruebas por debajo, considere la siguiente implementación manual:
# Implementación manual del test de White
#1. Construir el dataframe que relaciona el cuadrado del error con
#los términos lineales de la predictora y el término cuadrático
#(por ser modelo lineal simple no hay productos cruzados)
df_squared_error <- data.frame(
squared_error = mod1$residuals ^ 2,
precio_dolar = hass_dolar$precio_dolar,
precio_dolar_2 = hass_dolar$precio_dolar^2
)
#2. Calcule un modelo lineal ajustado sobre el cuadrado del error regresado por el termino lineal y cuadrático del precio del dolar
modelo_white <- lm(data = df_squared_error,
squared_error ~ precio_dolar + precio_dolar_2)
#3. Calcule el estadístico de prueba como n*R^2 (con n=datos del dataframe de errores)
R2 <- summary(modelo_white)$`r.squared`
w_estadistico = nrow(df_squared_error)*R2
#4. Como el estadístico se distribuye ji-cuadrado con 2 grados de libertad (2 variables regresoras) se tiene que el p valor es:
p_valor = 1-pchisq(as.numeric(w_estadistico),2)
paste("BP = ",round(w_estadistico,4)," y p_valor=",round(p_valor,5))## [1] "BP = 7.8925 y p_valor= 0.01933"
En la salida impresa puede notar que los resultados son consistentes con lo generado por la prueba de heterocedasticidad arrojada con la función bptest.
Además, cuando se estudia heterocedasticidad en segundo grado se aprecia que el modelo tiene varianzas del error no constantes porque estas dependen del cuadrado del precio del dólar. Veamos un poco esto:
estimated_squared_error <- modelo_white$coefficients[1] + modelo_white$coefficients[2]*df_squared_error$precio_dolar +modelo_white$coefficients[3]*df_squared_error$precio_dolar_2
p<- ggplot(data = df_squared_error, aes(x=precio_dolar,
y=squared_error))+
geom_point()+
geom_line(aes(y=estimated_squared_error, color="red"))+
labs(title="Modelo white sobre cuadrados del error",
subtitle = paste("R²=",round(R2,4),
"W-Estadístico=",round(w_estadistico,4),
"p_valor=",round(p_valor,5)))+
scale_color_manual(values=c("red"="blue"),
labels=c("Error cuadrado estimado"),
name="Modelo de White")
print(p)En este caso \(W = nR² = 7.8925\) es el estadístico de la prueba con \(k=2\) grados de libertad asociados a el número de covariables. Ahora bien \(\text{p valor} = 0.01933\), obtenido a través de calcular \(p(\chi^2_2>=7.8925)\), indicando esta vez que hay evidencia estadística para asumir heterocedasticidad de orden cuadrático. Podrá notar que un patrón parabólico no muy fuerte se observa cuando dibujamos el modelo sobre la nube de puntos del error cuadrado vs el precio del dólar.
Vamos a intentar resolver el problema de heterocedasticidad con una transformación de la variable respuesta que tienda a aplanar la curva de errores cuadrados responsable de la falla del supuesto. Dos posibles transformaciones son:
- Extraer raíz cuadrada de la respuesta.
- Extraer logaritmo natural de la respuesta.
La segunda transformación es mucho más severa. Probaremos con esta.
# Datos para nuevo modelo con y modificada a través de logaritmo natural
data_mod <- data.frame(
y_mod = log(hass_dolar$precio_aguacate_kg),
precio_dolar = hass_dolar$precio_dolar,
precio_dolar_2 = hass_dolar$precio_dolar^2
)
# Nuevo modelo lineal
mod2<-lm(y_mod~precio_dolar, data=data_mod)
residuales_mod<-mod2$residuals
#Guardamos los cuadrados del nuevo error
data_mod$squared_error <- residuales_mod^2
#Aplicamos prueba de white sobre el modelo con la y modificada
white_test_mod<-bptest(mod2,
varformula = ~ (precio_dolar + I(precio_dolar^2)),
data=data_mod)
print(white_test_mod)##
## studentized Breusch-Pagan test
##
## data: mod2
## BP = 5.2735, df = 2, p-value = 0.07159
w_estadistico_mod <- white_test_mod$statistic
p_value_mod <- white_test_mod$p.value
#Obtenemos los coeficientes para dibujar curva de errores cuadrados
#estimados
manual_white_test_mod <-lm(data = data_mod,
squared_error ~ precio_dolar + precio_dolar_2)
coefs_mod <- manual_white_test_mod$coefficients
R2_mod <- summary(manual_white_test_mod)$`r.squared`
#Obtenemos los nuevos errores cuadrados estimados
data_mod$estimated_squared_error <- (coefs_mod[1] +
coefs_mod[2]*data_mod$precio_dolar +
coefs_mod[3]*data_mod$precio_dolar^2)
#Comprobamos que no se nos deteriore la normalidad
Lilliefors_test_mod <- lillie.test(mod2$residuals)
print(Lilliefors_test_mod)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: mod2$residuals
## D = 0.050708, p-value = 0.544
#Comprobamos si milagrosamente resolvimos la no independencia del modelo
#original
resultado_ljung_box_mod <- Box.test(mod2$residuals, lag = i,
type = "Ljung-Box")
print(resultado_ljung_box_mod)##
## Box-Ljung test
##
## data: mod2$residuals
## X-squared = 133.25, df = 12, p-value < 0.00000000000000022
#Graficamos de nuevo el modelo de white sobre los cuadrados del error
#para el modelo con y modificada
p<- ggplot(data = data_mod, aes(x=precio_dolar,
y=squared_error))+
geom_point()+
geom_line(aes(y=estimated_squared_error, color="red"))+
labs(title="Modelo white sobre cuadrados del error",
subtitle = paste("R²=",round(R2_mod,4),
"W-Estadístico=",round(w_estadistico_mod,4),
"p_valor=",round(p_value_mod,5)))+
scale_color_manual(values=c("red"="blue"),
labels=c("Error cuadrado estimado"),
name="Modelo de White")
print(p)Las conclusiones de todo esto son:
- Nuestro modelo es inservible, porque la independencia sigue sin cumplirse, supuesto fundamental para poder usar el modelo con la confianza de no estar subestimando la varianza del error.
- La prueba de homocedasticidad pasa por un estrecho margen pese a las severas transformaciones.
- No tenemos problemas con normalidad.
Hasta acá las discusiones del modelo lineal simple con este dataset de ejemplo. Para nuestro capítulo siguiente relacionado con pruebas de hipótesis para los betas del modelo, tendremos que usar un modelo que cumpla con todos los supuestos. Para este efecto aprovecharemos para simular datos bajo un modelo teórico perfecto. Esta herramienta de simulación permitirá tomar control del cumplimiento de los supuestos para investigar tranquilamente el comportamiento estadístico de los betas estimados de un modelo cuando todos las condiciones para su uso están garantizadas.
Ejercicios propuestos:
Calcular un modelo de regresión para los cambios porcentuales del precio del aguacate de un mes a otro vs los cambios porcentuales del precio del dólar de un mes a otro. ¿En qué difiere este enfoque del originalmente seguido? ¿Cuál cree que sea la ventaja de generar un modelo de esta naturaleza?
Evalue la normalidad del error para este modelo usando la prueba de Lilliefors
Evalúe en un gráfico la consistencia de la prueba de Lilieforse (Pista: Use el código ya suministrado para calcular la máxima distancia entre la curva normal teórica y la curva empírica)
Pruebe la homocedasticidad del nuevo modelo usando tanto gráficos como pruebas formales, y dibuje la curva de errores cuadrados estimados basada en el modelo de Breusch-Pagan, es decir, sin incluir el término cuadrático. ¿Es esta prueba suficientemente concluyente?
Pruebe la independencia de los errores del modelo, comparéla contra los resultados obtenidos con el modelo que usaba los precios en valores absolutos, es decir, el modelo discutido en este documento.
Extraiga conclusiones generales de la adecuación del nuevo modelo a todos los supuestos, ¿Cuál de los dos modelos sería más útil para hacer predicciones, el modelo elaborado por los profesores o el que ustedes elaboraron? Felicítese mucho si el suyo es mejor. :)
Variable aleatoria: Una variable aleatoria X es una función de \[ f : \Omega \rightarrow \mathbb{R}\] con \(\Omega\) el conjunto de resultados posibles de un experimento aleatorio. En otras palabras un mapeo desde el espacio de descripciones muestrales a los números reales con el fin de asociar a los resultados del experimento un valor numérico con el cual habilitar operaciones matemáticas.
[Proceso estocástico]{style = “color: red;”}: Sea \(\mathcal{T}\) un conjunto arbitrario. Un proceso estocástico es una familia \(\{X(t), \text{ t} \in \tau \}\), tal que, para cada \(t \in \mathcal{T}, X(t)\) es una variable aleatoria. Es decir, es una secuencia de variables aleatorias ordenadas por un índice, típicamente tomado del conjunto de los números enteros, o del conjunto de los números reales.
[Proceso estocástico en tiempo discreto]{style = “color: red;}: Es un proceso estocástico con un índice temporal discreto. Significando que \(\tau \in \mathbb{Z}\)
[Especificación de un proceso estocástico]{style = “color: red;”}: Sea \(t_1, t_2, t_3,...,t_n\) elementos cualquiera de \(\mathcal{T}\) y consideremos
\[F(x_1, x_2, x_3, ..., x_n; t_1, t_2, t_3, ..., t_n)=P\{X(t_1)<=x_1, ..., X(t_n)<=x_n\}\] [Definición de un proceso estocástico estacionario]{style = “color: red;”}: Un proceso estocástico \(\{X(t), \text{ t} \in \mathcal{T} \}\) se dice estrictamente estacionario si todas las distribuciones finito dimensionales tal como se definen en la ecuación anterior permanecen siendo las mismas sobre traslaciones en el tiempo, es decir:
\[F(x_1, x_2, x_3, ..., x_n; t_1+\tau, t_2+\tau, t_3+\tau, ..., t_n+\tau)=F(x_1, x_2, x_3, ..., x_n; t_1, t_2, t_3, ..., t_n)\] Algunas consecuencias de lo anterior son:
- \(E\{X(t)\} = \mu(t) = \mu\)
- \(Var\{X(t)\} = \sigma²(t) = \sigma²\)
- \(\gamma(t_1, t_2) = \gamma(t_1-t_2, 0) = Cov\{X(t_1-t_2), X(0)\}\) o equivalentemente \(\gamma{\tau} = Cov\{X(t), X(t+\tau)\} = Cov\{X(0), X(\tau)\}\) para t, \(\tau \in \mathcal{T}\)
Reflexiones en lenguaje informal:
Un proceso estocástico quedará suficientemente definido si definimos la distribución de probabilidad n-dimensional conjunta para cada \(n>=1\), es decir, tomando \(n\) variables del proceso estocástico a la vez, sería necesario especificar la relación estadística de esas n variables a la vez para cada eleción de n. En nuestro contexto implica especificar cómo se espera que el precio de hoy se relacione con el precio de ayer, o cómo se espera que el precio de hoy se relacione con el de ayer y además con el de antier. Y así sucesivamente para cualquier grupo finito de precios que decida tomar a consideración.
Un proceso estocástico estacionario se puede interpretar como aquel cuyas propiedades estadísticas permanencen constantes a lo largo del tiempo y sólo cambiaran en función del grupo de \(n\) variables que decida obervar conjuntamente, pero no de la posición de esas n variables en el tiempo. Por ejemplo, suponiendo que n=1, esto implicaría que las propiedades estadísticas de cada variable tomada individualmente no se van alteradas conforme el tiempo avance, en particular cada variable será estable en media y en varianza. Lo que significa que puedo reunir las observaciones individuales en cada momento en el tiempo para estimar la media constante y la varianza de toda la secuencia, en la media en que en cada punto temporal estoy asumiendo que la media y la varianza es igual. Esto siempre se cumplirá? Por supuesto que no, es una suposición que imponemos a los datos para poder hacer inferencia con un único set de valores históricos observados.
Si la condición de estacionariedad estricta falla no hay nada que se pueda hacer en el estado actual del avance de la ciencia estadística.🤔 La única alternativa es relajar un poco las condiciones para exigencia de estacionariedad, cómo por ejemplo exigir estacionariedad débil consistente en pedir sólo homogeneidad en media, en varianza y en covarianza de segundo orden a la serie de tiempo bajo estudio. Exigir estacionariedad estricta es algo muy complicado de garantizar en la naturaleza.
Además es más fácil diseñar pruebas de estacionariedad débil.
Apliquemos una para la serie de tiempo de residuales
##
## Augmented Dickey-Fuller Test
##
## data: mod1$residuals
## Dickey-Fuller = -3.6959, Lag order = 5, p-value = 0.02729
## alternative hypothesis: stationary
# Prueba de White
resultado_white <- bptest(mod1)
# Muestra el resultado de la prueba
print(resultado_white)##
## studentized Breusch-Pagan test
##
## data: mod1
## BP = 1.741, df = 1, p-value = 0.187
En este caso debemos interpretar valores p<0.05 de la prueba de Dickey como un apoyo a la hipótesis alternativa de estacionariedad. Es decir, podemos considerar los residuales del modelo como estacionarios, porque un modelo autoregresivo con 5 rezagos no logró encontrar evidencia estadística de autocorrelaciones del error conducentes a no estacionariedad.
Por otro lado, podemos interpretar valores p<0.05 de la prueba de White como un apoyo a la hipótesis de homocedasticidad, es decir, podemos considerar los residuales del modelo como constantes en varianza a lo largo del tiempo.
Varianza de los residuales del modelo
Obteniendo los coeficientes.
## (Intercept) precio_dolar
## 1436.6790256 0.9449775
Valores ajustados - Fitted values.
## 1 2 3 4 5 6 7 8
## 3130.297 3123.597 3143.631 3139.479 3143.000 3157.035 3131.911 3109.503
## 9 10 11 12 13 14 15 16
## 3129.884 3148.096 3165.846 3183.871 3242.079 3230.649 3237.020 3248.871
## 17 18 19 20
## 3218.826 3253.241 3264.442 3291.015
Residuales
## 1 2 3 4 5 6 7
## 340.70260 335.90285 399.86948 652.52072 518.25021 423.46463 -185.91095
## 8 9 10 11 12 13 14
## -333.25250 -398.13417 164.70412 -127.34649 -238.37054 -245.07936 287.60071
## 15 16 17 18 19 20
## 581.57961 273.37874 246.42423 -73.04054 -439.44215 -171.26527
# Los valores ajustados y los residuales también se pueden recuperar usando las funciones fitted( ) y residuals( ). Consulte la ayuda de estas funciones para conocer otros detalles.Calcular la varianza de los residuales.
Porcentaje de variabilidad no explicada por el modelo.
Diagrama de dispersión con los puntos originales
Creación de la columna de predicción
En este chunk trabajamos con el modo de escritura Tidyverse, aunque se muestra cómo sería con R base.
#hass_dolar$predicciones <- predict(mod1)
hass_dolar <- hass_dolar %>%
mutate(predicciones = predict(mod1))hass_dolar %>%
ggplot(aes(x = precio_dolar, y = precio_aguacate_kg)) +
geom_smooth(method = "lm", se = FALSE, color="lightblue") +
geom_segment(aes(xend=precio_dolar, yend=predicciones),
col='red', lty='dashed') +
geom_point() +
geom_point(aes(y=predicciones), col='red') +
theme_light()REVISIÓN BÁSICA DE CONCEPTOS - Simulación.
Vamos a simular un modelo de regresión cuya especificación funcional es la siguiente:
\[ y = \beta_0 + \beta_1 * x + e \]
Con \(e \text{~} N(0,\sigma^2)\)
Como se puede notar, todo modelo teorico se compone de dos términos:
- El valor E(y|x)= E_y_x (valor esperado de y dado x)
- El componente aleatorio también llamado error.
Aunque en el modelo anterior hemos elegido una estructura lineal para representar E(y|x) en realidad podemos elegir cualquier otra estructura alternativa, siempre que respetemos la linealidad en los betas.
bo = 2
b1_1 = 3
b1_2 = 0.5
x <- seq(1,10)
# Estructura para el valor esperado de y dado x ()
# Esta estructura admitiría otras formas funcionales
# El valor esperado es deterministico.
E_y_x_1 <- bo + b1_1 * x
E_y_x_2 <- bo + b1_1 * x + b1_2 * x^2
# el valor observado o valor real, es aleatorio.
y_obs_1 <- E_y_x_1 + rnorm(10, 0, 4)
y_obs_2 <- E_y_x_2 + rnorm(10, 0, 4)
datos_simulados <- data.frame(
x = x,
x_2 = x^2,
E_y_x_1 = E_y_x_1,
E_y_x_2 = E_y_x_2,
y_obs_1 = y_obs_1,
y_obs_2 = y_obs_2,
y_pred_1 = predict(lm(y_obs_1 ~ x)),
y_pred_2 = predict(lm(y_obs_2 ~ x+ x^2)))Dibujamos un gráfico de dispersión
datos_simulados %>%
ggplot(aes(x = x, y = y_obs_1)) +
geom_smooth(method = "lm", formula = y~x, se = FALSE, color="lightblue") +
geom_point(col = "green")+
geom_point(aes(y=y_obs_2), col='red')+
geom_smooth(method="lm", se= FALSE ,formula=y_obs_2~poly(x, 2),color = "blue")+
ylab("Y")DISCUTIENDO SOBRE LA LINEALIDAD DE LOS BETAS
En esta sección vamos a aclarar el asunto de que los BETAS sean lineales.
La linealidad en este contexto hace referencia a que no tengan potencias y que los coeficientes sean independendientes entre ellos.
Supongamos un modelo con coeficientes repetidos.
\[E(y|x) = \beta_0 + \alpha * x + \alpha * x^2 \]
Debemos factorizarlos para evitar estimarlos por separado y producir inconsistencias,esto debido a que ambos alpha (coeficientes) son dependientes.Es decir, reescribimos:
\[E(y|x) = \beta_0 + \alpha (x + x^2) \]
Donde el componente x + x² será el vector de valores de la variable predictora, que puede pensarse como una variable nueva:
\[ \tilde{x} = x + x^2\]
Y por lo tanto pensar al modelo como:
\[ E(y|x) = \beta_0 + \alpha * \tilde{x}\]
Esta transformación garantiza la creación de un nuevo modelo lineal respecto a \(\beta_0\) y \(\alpha\).
Es importante anotar que el dataframe que alimentará el modelo deberá tener computada la variable auxiliar \(\tilde{x}\), la cual debe ser pasada a la función encargada de estimar parámetros.
Esto corresponde a crear una nueva variable con los cálculos mencionados \(\tilde{x} = x + x^2\). Sin embargo, para graficar es conveniente usar la variable original.
# Se define una semilla para generar los aleatorios.
set.seed(77)
alpha <- 0.7
E_y_x_3 = bo + alpha * (x + x^2)
y_obs_3 = E_y_x_3 + rnorm(10, 0, 4)
x_auxiliar = x + x^2
datos_simulados_2 <- data.frame(
x = x,
x_2 = x^2,
x_auxiliar = x_auxiliar,
E_y_x_3 = E_y_x_3,
y_obs_3 = y_obs_3,
y_pred_3_1 = predict(lm(y_obs_3 ~ poly(x, 2))),
y_pred_3_2 = predict(lm(y_obs_3 ~ x_auxiliar )))
modelo1a <- lm(y_obs_3 ~ poly(x, 2))
modelo1b <- lm(y_obs_3 ~ poly(x, 2, raw = TRUE))
modelo2 <- lm(y_obs_3 ~ x_auxiliar )Notese que los coeficientes arrojados por el modelo1a al invocar la función poly sin parámetro de raw=TRUE, no son atribuibles o no están asociados a los términos \(1,x,x^2\), es decir, No es válido escribir:
\[E(y|x)= 34.41 * 1 + 76.82 * x + 16.74 * x^2\]
Esta ecuación no sería válida porque al remplazar para \(x=5\) se obtiene \(E(y|5)=\) 836.91, mientras que al remplazar en el modelo con el parámetro raw = TRUE, se obtendría: \(E(y|5)=\) 24.35.
Observe que para el valor \(x=5\) los valores observados de \(y\) están cercanos a 24 por lo que una predicción arrojando un valor cercano a 836, está claramente desfasada. Observe el gráfico más abajo para comprender lo explicado.
Esto se debe a que los coeficientes anclados al modelo_1a están diseñados para ser usados sobre la siguiente base de polinomios ortogonales:
- \(P_1(x)= a\)
- \(P_2(x)=m*x - c\)
- \(P_3(x)= k_1 + k_2 *x + k_3 * x ^2\)
Donde \(a,m,c,k_1,k_2,k_3\) son parámetros estimados que lamentablemente no están siendo arrojados por el modelo.
A continuación se muestra el resumen de cada uno de los modelos, utilizando la función summary().
##
## Call:
## lm(formula = y_obs_3 ~ poly(x, 2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4051 -2.0000 0.2405 2.8394 3.5594
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.407 1.143 30.108 0.0000000115 ***
## poly(x, 2)1 76.820 3.614 21.257 0.0000001284 ***
## poly(x, 2)2 16.736 3.614 4.631 0.00239 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.614 on 7 degrees of freedom
## Multiple R-squared: 0.9854, Adjusted R-squared: 0.9813
## F-statistic: 236.7 on 2 and 7 DF, p-value: 0.0000003736
##
## Call:
## lm(formula = y_obs_3 ~ poly(x, 2, raw = TRUE))
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4051 -2.0000 0.2405 2.8394 3.5594
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.9133 4.2503 0.921 0.38783
## poly(x, 2, raw = TRUE)1 0.4458 1.7751 0.251 0.80890
## poly(x, 2, raw = TRUE)2 0.7283 0.1573 4.631 0.00239 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.614 on 7 degrees of freedom
## Multiple R-squared: 0.9854, Adjusted R-squared: 0.9813
## F-statistic: 236.7 on 2 and 7 DF, p-value: 0.0000003736
##
## Call:
## lm(formula = y_obs_3 ~ x_auxiliar)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5584 -2.0443 0.1458 2.8929 3.7758
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.35135 1.71306 1.956 0.0861 .
## x_auxiliar 0.70580 0.03039 23.222 0.0000000126 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.386 on 8 degrees of freedom
## Multiple R-squared: 0.9854, Adjusted R-squared: 0.9836
## F-statistic: 539.3 on 1 and 8 DF, p-value: 0.00000001256
datos_simulados_2 %>%
ggplot(aes(x = x , y = y_obs_3))+
geom_point(col = 'green')+
geom_point(aes(y = y_pred_3_1), col = 'blue')+
geom_point(aes(y = y_pred_3_2), col = 'red')ESTIMACIÓN DE PARÁMETROS
Por máxima verosimilitud.
Por mínimos cuadrados.
n <- 1000
vector_x <- c(1: n )
sigma <- 4
beta_0 <- 2
beta_1 <- 0.5
# Esperada de y dado x
E_y_x <- beta_0 + (beta_1 * vector_x )
print(E_y_x)## [1] 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0
## [13] 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0
## [25] 14.5 15.0 15.5 16.0 16.5 17.0 17.5 18.0 18.5 19.0 19.5 20.0
## [37] 20.5 21.0 21.5 22.0 22.5 23.0 23.5 24.0 24.5 25.0 25.5 26.0
## [49] 26.5 27.0 27.5 28.0 28.5 29.0 29.5 30.0 30.5 31.0 31.5 32.0
## [61] 32.5 33.0 33.5 34.0 34.5 35.0 35.5 36.0 36.5 37.0 37.5 38.0
## [73] 38.5 39.0 39.5 40.0 40.5 41.0 41.5 42.0 42.5 43.0 43.5 44.0
## [85] 44.5 45.0 45.5 46.0 46.5 47.0 47.5 48.0 48.5 49.0 49.5 50.0
## [97] 50.5 51.0 51.5 52.0 52.5 53.0 53.5 54.0 54.5 55.0 55.5 56.0
## [109] 56.5 57.0 57.5 58.0 58.5 59.0 59.5 60.0 60.5 61.0 61.5 62.0
## [121] 62.5 63.0 63.5 64.0 64.5 65.0 65.5 66.0 66.5 67.0 67.5 68.0
## [133] 68.5 69.0 69.5 70.0 70.5 71.0 71.5 72.0 72.5 73.0 73.5 74.0
## [145] 74.5 75.0 75.5 76.0 76.5 77.0 77.5 78.0 78.5 79.0 79.5 80.0
## [157] 80.5 81.0 81.5 82.0 82.5 83.0 83.5 84.0 84.5 85.0 85.5 86.0
## [169] 86.5 87.0 87.5 88.0 88.5 89.0 89.5 90.0 90.5 91.0 91.5 92.0
## [181] 92.5 93.0 93.5 94.0 94.5 95.0 95.5 96.0 96.5 97.0 97.5 98.0
## [193] 98.5 99.0 99.5 100.0 100.5 101.0 101.5 102.0 102.5 103.0 103.5 104.0
## [205] 104.5 105.0 105.5 106.0 106.5 107.0 107.5 108.0 108.5 109.0 109.5 110.0
## [217] 110.5 111.0 111.5 112.0 112.5 113.0 113.5 114.0 114.5 115.0 115.5 116.0
## [229] 116.5 117.0 117.5 118.0 118.5 119.0 119.5 120.0 120.5 121.0 121.5 122.0
## [241] 122.5 123.0 123.5 124.0 124.5 125.0 125.5 126.0 126.5 127.0 127.5 128.0
## [253] 128.5 129.0 129.5 130.0 130.5 131.0 131.5 132.0 132.5 133.0 133.5 134.0
## [265] 134.5 135.0 135.5 136.0 136.5 137.0 137.5 138.0 138.5 139.0 139.5 140.0
## [277] 140.5 141.0 141.5 142.0 142.5 143.0 143.5 144.0 144.5 145.0 145.5 146.0
## [289] 146.5 147.0 147.5 148.0 148.5 149.0 149.5 150.0 150.5 151.0 151.5 152.0
## [301] 152.5 153.0 153.5 154.0 154.5 155.0 155.5 156.0 156.5 157.0 157.5 158.0
## [313] 158.5 159.0 159.5 160.0 160.5 161.0 161.5 162.0 162.5 163.0 163.5 164.0
## [325] 164.5 165.0 165.5 166.0 166.5 167.0 167.5 168.0 168.5 169.0 169.5 170.0
## [337] 170.5 171.0 171.5 172.0 172.5 173.0 173.5 174.0 174.5 175.0 175.5 176.0
## [349] 176.5 177.0 177.5 178.0 178.5 179.0 179.5 180.0 180.5 181.0 181.5 182.0
## [361] 182.5 183.0 183.5 184.0 184.5 185.0 185.5 186.0 186.5 187.0 187.5 188.0
## [373] 188.5 189.0 189.5 190.0 190.5 191.0 191.5 192.0 192.5 193.0 193.5 194.0
## [385] 194.5 195.0 195.5 196.0 196.5 197.0 197.5 198.0 198.5 199.0 199.5 200.0
## [397] 200.5 201.0 201.5 202.0 202.5 203.0 203.5 204.0 204.5 205.0 205.5 206.0
## [409] 206.5 207.0 207.5 208.0 208.5 209.0 209.5 210.0 210.5 211.0 211.5 212.0
## [421] 212.5 213.0 213.5 214.0 214.5 215.0 215.5 216.0 216.5 217.0 217.5 218.0
## [433] 218.5 219.0 219.5 220.0 220.5 221.0 221.5 222.0 222.5 223.0 223.5 224.0
## [445] 224.5 225.0 225.5 226.0 226.5 227.0 227.5 228.0 228.5 229.0 229.5 230.0
## [457] 230.5 231.0 231.5 232.0 232.5 233.0 233.5 234.0 234.5 235.0 235.5 236.0
## [469] 236.5 237.0 237.5 238.0 238.5 239.0 239.5 240.0 240.5 241.0 241.5 242.0
## [481] 242.5 243.0 243.5 244.0 244.5 245.0 245.5 246.0 246.5 247.0 247.5 248.0
## [493] 248.5 249.0 249.5 250.0 250.5 251.0 251.5 252.0 252.5 253.0 253.5 254.0
## [505] 254.5 255.0 255.5 256.0 256.5 257.0 257.5 258.0 258.5 259.0 259.5 260.0
## [517] 260.5 261.0 261.5 262.0 262.5 263.0 263.5 264.0 264.5 265.0 265.5 266.0
## [529] 266.5 267.0 267.5 268.0 268.5 269.0 269.5 270.0 270.5 271.0 271.5 272.0
## [541] 272.5 273.0 273.5 274.0 274.5 275.0 275.5 276.0 276.5 277.0 277.5 278.0
## [553] 278.5 279.0 279.5 280.0 280.5 281.0 281.5 282.0 282.5 283.0 283.5 284.0
## [565] 284.5 285.0 285.5 286.0 286.5 287.0 287.5 288.0 288.5 289.0 289.5 290.0
## [577] 290.5 291.0 291.5 292.0 292.5 293.0 293.5 294.0 294.5 295.0 295.5 296.0
## [589] 296.5 297.0 297.5 298.0 298.5 299.0 299.5 300.0 300.5 301.0 301.5 302.0
## [601] 302.5 303.0 303.5 304.0 304.5 305.0 305.5 306.0 306.5 307.0 307.5 308.0
## [613] 308.5 309.0 309.5 310.0 310.5 311.0 311.5 312.0 312.5 313.0 313.5 314.0
## [625] 314.5 315.0 315.5 316.0 316.5 317.0 317.5 318.0 318.5 319.0 319.5 320.0
## [637] 320.5 321.0 321.5 322.0 322.5 323.0 323.5 324.0 324.5 325.0 325.5 326.0
## [649] 326.5 327.0 327.5 328.0 328.5 329.0 329.5 330.0 330.5 331.0 331.5 332.0
## [661] 332.5 333.0 333.5 334.0 334.5 335.0 335.5 336.0 336.5 337.0 337.5 338.0
## [673] 338.5 339.0 339.5 340.0 340.5 341.0 341.5 342.0 342.5 343.0 343.5 344.0
## [685] 344.5 345.0 345.5 346.0 346.5 347.0 347.5 348.0 348.5 349.0 349.5 350.0
## [697] 350.5 351.0 351.5 352.0 352.5 353.0 353.5 354.0 354.5 355.0 355.5 356.0
## [709] 356.5 357.0 357.5 358.0 358.5 359.0 359.5 360.0 360.5 361.0 361.5 362.0
## [721] 362.5 363.0 363.5 364.0 364.5 365.0 365.5 366.0 366.5 367.0 367.5 368.0
## [733] 368.5 369.0 369.5 370.0 370.5 371.0 371.5 372.0 372.5 373.0 373.5 374.0
## [745] 374.5 375.0 375.5 376.0 376.5 377.0 377.5 378.0 378.5 379.0 379.5 380.0
## [757] 380.5 381.0 381.5 382.0 382.5 383.0 383.5 384.0 384.5 385.0 385.5 386.0
## [769] 386.5 387.0 387.5 388.0 388.5 389.0 389.5 390.0 390.5 391.0 391.5 392.0
## [781] 392.5 393.0 393.5 394.0 394.5 395.0 395.5 396.0 396.5 397.0 397.5 398.0
## [793] 398.5 399.0 399.5 400.0 400.5 401.0 401.5 402.0 402.5 403.0 403.5 404.0
## [805] 404.5 405.0 405.5 406.0 406.5 407.0 407.5 408.0 408.5 409.0 409.5 410.0
## [817] 410.5 411.0 411.5 412.0 412.5 413.0 413.5 414.0 414.5 415.0 415.5 416.0
## [829] 416.5 417.0 417.5 418.0 418.5 419.0 419.5 420.0 420.5 421.0 421.5 422.0
## [841] 422.5 423.0 423.5 424.0 424.5 425.0 425.5 426.0 426.5 427.0 427.5 428.0
## [853] 428.5 429.0 429.5 430.0 430.5 431.0 431.5 432.0 432.5 433.0 433.5 434.0
## [865] 434.5 435.0 435.5 436.0 436.5 437.0 437.5 438.0 438.5 439.0 439.5 440.0
## [877] 440.5 441.0 441.5 442.0 442.5 443.0 443.5 444.0 444.5 445.0 445.5 446.0
## [889] 446.5 447.0 447.5 448.0 448.5 449.0 449.5 450.0 450.5 451.0 451.5 452.0
## [901] 452.5 453.0 453.5 454.0 454.5 455.0 455.5 456.0 456.5 457.0 457.5 458.0
## [913] 458.5 459.0 459.5 460.0 460.5 461.0 461.5 462.0 462.5 463.0 463.5 464.0
## [925] 464.5 465.0 465.5 466.0 466.5 467.0 467.5 468.0 468.5 469.0 469.5 470.0
## [937] 470.5 471.0 471.5 472.0 472.5 473.0 473.5 474.0 474.5 475.0 475.5 476.0
## [949] 476.5 477.0 477.5 478.0 478.5 479.0 479.5 480.0 480.5 481.0 481.5 482.0
## [961] 482.5 483.0 483.5 484.0 484.5 485.0 485.5 486.0 486.5 487.0 487.5 488.0
## [973] 488.5 489.0 489.5 490.0 490.5 491.0 491.5 492.0 492.5 493.0 493.5 494.0
## [985] 494.5 495.0 495.5 496.0 496.5 497.0 497.5 498.0 498.5 499.0 499.5 500.0
## [997] 500.5 501.0 501.5 502.0
error <- rnorm(n = n, mean = 0, sd = sigma)
y_observado <- beta_0 + (beta_1 * vector_x ) + error
print(y_observado)## [1] -9.265541 2.028646 2.937660 3.869384 5.619264 7.360592
## [7] 9.597181 14.429291 7.118537 10.652279 6.483102 14.077727
## [13] 15.625079 5.484622 3.383127 10.543646 7.664401 5.361710
## [19] 18.823746 17.160776 3.050258 10.796101 12.281348 11.000155
## [25] 15.074887 12.802274 16.139696 15.649304 16.824517 20.596477
## [31] 17.512477 15.875106 15.660199 17.834672 23.038678 19.382678
## [37] 16.682756 23.669010 23.053954 24.196648 17.280075 26.547745
## [43] 32.843871 26.012656 15.428101 16.873824 24.947705 22.189917
## [49] 32.810645 22.008560 29.094688 28.106801 24.494048 29.035585
## [55] 28.843687 31.690343 28.804476 29.121890 32.188607 39.387922
## [61] 25.946999 31.670173 35.291761 35.089666 33.482570 44.393269
## [67] 34.647942 36.202546 27.516022 40.678747 40.505934 37.733686
## [73] 36.485893 38.358421 39.127831 41.895726 41.022360 43.869311
## [79] 41.487409 33.293631 46.007014 46.340168 43.307650 30.361798
## [85] 38.446858 46.009403 43.606182 49.908762 56.240348 46.670516
## [91] 44.916267 41.514288 49.015419 46.212810 42.992990 48.778560
## [97] 51.668788 55.489479 48.413265 53.911428 55.407637 53.495620
## [103] 49.383730 51.264664 53.980522 55.049370 52.070863 61.650524
## [109] 56.008867 61.766359 60.503476 64.118084 53.906824 65.393542
## [115] 52.102992 60.585501 64.699350 63.053040 62.412397 65.718989
## [121] 56.698966 62.316500 58.716476 61.288476 64.163728 66.896376
## [127] 68.362137 67.107193 73.120807 61.198593 64.914115 61.154078
## [133] 68.453434 73.335860 69.347666 65.677687 71.764833 64.271905
## [139] 68.477214 74.473414 75.625322 73.046914 74.782204 75.483451
## [145] 78.939615 81.572851 74.979472 66.748252 74.872171 79.686641
## [151] 76.264237 77.114668 76.972465 75.874119 74.614837 79.883959
## [157] 77.197678 77.969892 79.334054 76.600718 80.710794 85.912116
## [163] 82.575927 80.249849 81.895177 80.493869 84.603712 81.021345
## [169] 77.383178 82.125694 86.757934 92.194819 79.615810 88.402900
## [175] 90.099991 95.780725 95.664267 91.660829 92.307321 91.199442
## [181] 94.751780 90.627674 95.740870 92.186361 97.273126 94.751578
## [187] 90.615784 98.208120 95.274842 99.066993 93.981242 96.985319
## [193] 100.018279 96.917646 96.605909 105.875634 102.882755 105.208926
## [199] 104.852755 102.995486 98.944024 100.582438 105.375775 109.599393
## [205] 102.171532 109.387365 104.906178 109.087417 103.100220 103.003452
## [211] 112.676571 104.410939 112.737950 109.335340 107.930901 105.857007
## [217] 105.491265 107.616343 110.896275 108.490330 108.247700 118.422837
## [223] 116.485316 111.456689 119.010276 109.286637 114.700014 113.242019
## [229] 114.994203 120.846571 115.551315 116.582148 120.711077 120.278320
## [235] 122.503374 113.266525 118.014070 116.035079 121.798163 118.172809
## [241] 127.292661 128.772185 127.245496 122.553616 122.868983 123.767254
## [247] 123.603193 124.068996 129.793829 129.567219 129.013043 129.647393
## [253] 129.622391 127.337691 134.011768 128.653053 137.296293 131.022219
## [259] 130.757112 128.646890 131.515937 127.278725 141.599687 137.645591
## [265] 141.968057 137.672219 133.688976 138.777617 134.443015 138.013143
## [271] 134.369128 133.671307 128.625288 140.630644 145.346955 144.802831
## [277] 138.823189 143.667246 144.078795 146.514479 143.870761 143.596282
## [283] 143.613048 146.037859 150.991206 150.534877 150.592148 149.807308
## [289] 148.593098 149.533588 144.753246 155.056930 151.141460 153.192276
## [295] 149.782141 145.637842 151.104023 154.776958 149.998783 151.083057
## [301] 155.671462 153.597301 152.461846 152.916740 151.863383 156.845020
## [307] 154.530705 162.427797 155.602518 161.750836 153.690598 161.469722
## [313] 163.010299 153.044913 159.538016 166.086049 158.597632 163.229112
## [319] 158.824643 167.615507 163.395024 167.977319 166.728506 170.109100
## [325] 159.285369 167.169306 160.694178 166.571654 164.073120 157.655112
## [331] 168.798608 163.909142 167.974416 172.856515 175.650578 176.555039
## [337] 167.631605 172.972786 162.076239 170.648110 178.404094 173.559135
## [343] 167.206168 173.048448 169.069481 175.549638 177.176409 170.474740
## [349] 175.798856 177.497981 171.247445 173.889131 173.400347 180.328334
## [355] 180.342073 176.476584 188.526120 181.619194 191.496570 181.262353
## [361] 179.893624 184.209107 185.332254 180.102931 176.275280 191.259255
## [367] 191.665894 185.236948 185.286964 186.766325 190.446317 184.653539
## [373] 201.359483 192.627786 189.034749 187.805433 189.639361 191.896128
## [379] 195.750745 187.720068 189.751628 188.643408 199.948917 192.505760
## [385] 191.718299 197.306856 197.712207 196.477588 195.991698 203.535596
## [391] 200.707341 196.576718 199.497179 201.204693 195.567848 199.331934
## [397] 207.241355 201.144083 202.815188 200.133860 196.282424 202.709178
## [403] 206.811808 206.152802 199.338514 208.203747 210.446042 210.714271
## [409] 213.304216 204.601236 208.201704 209.413339 214.802023 215.015860
## [415] 208.320534 208.956146 207.300618 215.920216 214.535939 208.922329
## [421] 214.204878 211.946548 213.509210 217.891529 217.815421 211.381764
## [427] 219.141138 214.824631 215.151008 218.092296 219.110325 216.326643
## [433] 220.660202 214.694054 227.341562 219.280536 223.373807 223.593935
## [439] 222.583633 223.747752 225.602931 224.693224 218.941212 229.890763
## [445] 217.544335 220.381504 230.224435 227.171083 227.023340 222.877020
## [451] 227.791909 227.403662 227.841227 226.364997 229.234458 225.752192
## [457] 227.598351 231.206864 230.960630 241.206867 233.941380 232.531573
## [463] 237.738696 239.908031 230.863586 238.311611 238.307978 238.798197
## [469] 234.331007 236.200720 237.463167 234.832624 239.585767 237.694200
## [475] 234.176699 234.139977 243.225557 243.056163 243.671284 241.546864
## [481] 244.544265 249.242251 248.957388 248.034055 245.673904 246.375226
## [487] 239.318235 244.312557 247.673404 245.201081 256.701859 245.913761
## [493] 242.644201 250.114914 248.624974 249.976235 251.013811 251.194173
## [499] 252.460027 244.396518 261.673313 258.818199 254.874468 249.011531
## [505] 253.925723 255.987432 251.605316 251.257305 260.143437 259.215913
## [511] 258.256589 260.647611 263.002572 259.742351 253.883285 259.526766
## [517] 259.652476 263.559505 260.076776 263.746374 258.182144 261.884827
## [523] 263.866577 266.771189 260.547230 263.594274 270.642094 272.201849
## [529] 267.419493 264.162520 265.438920 269.742620 264.685308 272.030936
## [535] 266.662131 269.108769 277.434810 272.760171 274.843552 273.864145
## [541] 282.742115 269.354760 275.599835 277.247630 281.571713 272.139505
## [547] 274.436391 277.355075 274.713388 283.432214 279.420251 286.114388
## [553] 282.348281 273.943627 280.473689 276.316526 285.247855 277.284674
## [559] 277.565818 283.746303 279.572771 285.156227 284.650010 279.473353
## [565] 288.668110 285.515644 285.986894 283.499818 289.623779 289.381117
## [571] 287.574424 289.557505 294.374395 279.904182 296.895110 293.219946
## [577] 288.596349 299.389704 290.404111 299.200482 294.524607 303.924934
## [583] 290.762378 289.362536 291.827077 295.611421 297.516090 293.703960
## [589] 295.411456 293.375499 295.589814 297.079774 297.640001 301.602918
## [595] 299.911704 295.834171 290.466146 300.997352 305.775667 302.230375
## [601] 304.671462 299.933488 311.105334 304.726298 308.549661 304.235342
## [607] 308.231526 307.847146 300.090793 314.325897 308.148647 303.713190
## [613] 307.551379 302.116400 314.127314 309.892431 309.381423 309.104443
## [619] 313.557965 311.586999 312.129004 308.506769 312.931923 312.953436
## [625] 316.331285 324.177663 318.382789 320.756939 319.920987 318.994367
## [631] 317.542702 314.234985 321.420195 325.704414 316.645132 316.044825
## [637] 323.675795 321.556617 315.955763 313.861739 326.795536 325.736815
## [643] 314.824302 328.099728 325.451974 325.131236 322.632937 319.990054
## [649] 323.996140 327.767415 328.114672 329.152268 330.817721 331.024415
## [655] 328.242358 328.278875 330.653693 329.507130 325.828999 326.192650
## [661] 337.087009 327.787032 337.225368 331.406881 332.791433 335.312026
## [667] 337.022039 336.536775 339.075566 338.160727 338.810382 339.692456
## [673] 347.123512 334.282048 339.347978 345.152070 337.340515 341.639172
## [679] 337.794137 333.896360 342.380822 348.802788 347.337464 348.013165
## [685] 344.282128 346.433430 347.921982 342.710095 350.147075 356.782976
## [691] 344.890280 349.608202 353.883495 343.115519 349.996676 351.305610
## [697] 352.609894 353.368425 352.224441 353.846822 352.805416 352.789138
## [703] 355.843507 350.428624 347.651100 359.808707 356.012044 360.301229
## [709] 356.181003 362.608294 358.108787 352.051227 353.298969 360.879583
## [715] 357.689598 360.013792 362.437046 352.775901 358.151468 367.812651
## [721] 369.767543 359.105722 359.862524 363.146111 367.289517 364.445608
## [727] 367.455366 361.254291 369.916434 370.974973 372.888598 365.720296
## [733] 372.749582 373.363895 375.130354 369.862537 367.847975 378.489652
## [739] 369.590009 368.871551 371.554602 369.904352 371.130991 370.818494
## [745] 372.861006 377.920785 375.651088 371.464377 365.456164 378.735999
## [751] 374.112540 376.222674 380.117902 375.737087 385.965530 382.047476
## [757] 376.267920 375.594777 377.486834 380.990322 385.345196 381.988560
## [763] 382.431276 378.610329 384.868176 393.130869 380.631748 391.371372
## [769] 385.548957 390.697413 390.585586 387.904856 390.087681 389.511307
## [775] 389.380056 387.856161 391.634503 394.300604 396.806255 379.855424
## [781] 390.471544 400.809673 396.711576 395.383267 389.233582 392.300110
## [787] 401.094097 401.381000 391.341719 392.424322 395.317025 396.663795
## [793] 395.351615 401.725386 397.213723 401.180504 399.272330 401.351139
## [799] 401.700503 400.114878 408.131808 400.057449 405.242381 407.597894
## [805] 402.652567 403.091423 406.806522 400.185329 403.455906 406.794075
## [811] 409.950616 402.533425 399.461950 407.674106 410.738648 406.507409
## [817] 407.676525 407.364338 408.216298 422.904292 413.673786 413.242213
## [823] 415.496754 408.188281 409.254191 415.178024 416.157409 417.336312
## [829] 414.591271 411.960069 417.008112 421.149164 418.481538 422.596943
## [835] 420.405217 420.579384 419.611500 420.794681 422.380343 416.232951
## [841] 420.503869 427.111943 421.115965 425.895258 425.958656 425.168995
## [847] 419.743363 420.193497 425.488467 419.938266 426.745067 428.146736
## [853] 427.314940 419.705481 435.009791 433.417581 430.169851 429.229599
## [859] 427.642260 427.291852 429.187161 432.963893 438.743576 438.066404
## [865] 431.010002 437.071490 432.797783 431.654688 435.738652 438.690006
## [871] 436.654107 446.120207 437.926523 435.403460 435.738268 443.334806
## [877] 444.917374 441.703131 445.001219 438.250156 438.160136 439.373964
## [883] 447.093405 438.355320 440.106577 446.414397 447.145599 448.760267
## [889] 444.726729 448.596649 439.658748 446.982974 453.676815 450.145542
## [895] 449.559567 457.948035 456.097105 443.728433 449.902008 449.645858
## [901] 449.559485 453.252741 451.738857 458.768271 459.618346 455.813839
## [907] 451.579433 448.664233 457.347313 453.468386 456.707244 460.500436
## [913] 454.882857 454.752996 458.145468 461.229505 457.616242 451.693008
## [919] 466.888626 459.676032 463.005717 467.287189 458.914782 462.929747
## [925] 469.151031 469.746548 463.342104 466.386190 458.349465 464.058421
## [931] 468.730267 466.582424 467.383085 466.329555 471.921350 473.918364
## [937] 472.374704 470.131541 474.473193 469.085841 471.435022 474.532248
## [943] 473.337380 471.497225 475.645770 475.016123 472.348478 477.556163
## [949] 481.692415 475.382491 479.442664 478.100581 487.983723 478.771865
## [955] 483.452833 478.191372 478.773329 486.572550 474.961427 483.150533
## [961] 489.676872 482.823038 487.571504 482.421791 484.032017 481.828466
## [967] 485.024457 488.403014 483.563388 485.432581 487.674154 484.832976
## [973] 489.882549 488.167992 489.861754 491.503086 483.244408 498.989761
## [979] 490.308198 489.267429 492.533820 489.513076 497.123950 496.822993
## [985] 495.545213 492.637233 489.568063 499.809256 493.638623 490.173591
## [991] 498.454394 502.333390 498.023531 504.804083 499.643065 503.938944
## [997] 503.038344 505.280834 499.543430 498.079409
Modelo lineal
modelo_y_predicho <- lm(y_observado ~ vector_x)
y_predicho <- predict(modelo_y_predicho)
print(y_predicho)## 1 2 3 4 5 6 7
## 2.650034 3.149898 3.649763 4.149627 4.649491 5.149356 5.649220
## 8 9 10 11 12 13 14
## 6.149084 6.648949 7.148813 7.648678 8.148542 8.648406 9.148271
## 15 16 17 18 19 20 21
## 9.648135 10.148000 10.647864 11.147728 11.647593 12.147457 12.647322
## 22 23 24 25 26 27 28
## 13.147186 13.647050 14.146915 14.646779 15.146643 15.646508 16.146372
## 29 30 31 32 33 34 35
## 16.646237 17.146101 17.645965 18.145830 18.645694 19.145559 19.645423
## 36 37 38 39 40 41 42
## 20.145287 20.645152 21.145016 21.644881 22.144745 22.644609 23.144474
## 43 44 45 46 47 48 49
## 23.644338 24.144202 24.644067 25.143931 25.643796 26.143660 26.643524
## 50 51 52 53 54 55 56
## 27.143389 27.643253 28.143118 28.642982 29.142846 29.642711 30.142575
## 57 58 59 60 61 62 63
## 30.642440 31.142304 31.642168 32.142033 32.641897 33.141761 33.641626
## 64 65 66 67 68 69 70
## 34.141490 34.641355 35.141219 35.641083 36.140948 36.640812 37.140677
## 71 72 73 74 75 76 77
## 37.640541 38.140405 38.640270 39.140134 39.639999 40.139863 40.639727
## 78 79 80 81 82 83 84
## 41.139592 41.639456 42.139320 42.639185 43.139049 43.638914 44.138778
## 85 86 87 88 89 90 91
## 44.638642 45.138507 45.638371 46.138236 46.638100 47.137964 47.637829
## 92 93 94 95 96 97 98
## 48.137693 48.637558 49.137422 49.637286 50.137151 50.637015 51.136879
## 99 100 101 102 103 104 105
## 51.636744 52.136608 52.636473 53.136337 53.636201 54.136066 54.635930
## 106 107 108 109 110 111 112
## 55.135795 55.635659 56.135523 56.635388 57.135252 57.635117 58.134981
## 113 114 115 116 117 118 119
## 58.634845 59.134710 59.634574 60.134438 60.634303 61.134167 61.634032
## 120 121 122 123 124 125 126
## 62.133896 62.633760 63.133625 63.633489 64.133354 64.633218 65.133082
## 127 128 129 130 131 132 133
## 65.632947 66.132811 66.632676 67.132540 67.632404 68.132269 68.632133
## 134 135 136 137 138 139 140
## 69.131997 69.631862 70.131726 70.631591 71.131455 71.631319 72.131184
## 141 142 143 144 145 146 147
## 72.631048 73.130913 73.630777 74.130641 74.630506 75.130370 75.630235
## 148 149 150 151 152 153 154
## 76.130099 76.629963 77.129828 77.629692 78.129556 78.629421 79.129285
## 155 156 157 158 159 160 161
## 79.629150 80.129014 80.628878 81.128743 81.628607 82.128472 82.628336
## 162 163 164 165 166 167 168
## 83.128200 83.628065 84.127929 84.627794 85.127658 85.627522 86.127387
## 169 170 171 172 173 174 175
## 86.627251 87.127115 87.626980 88.126844 88.626709 89.126573 89.626437
## 176 177 178 179 180 181 182
## 90.126302 90.626166 91.126031 91.625895 92.125759 92.625624 93.125488
## 183 184 185 186 187 188 189
## 93.625352 94.125217 94.625081 95.124946 95.624810 96.124674 96.624539
## 190 191 192 193 194 195 196
## 97.124403 97.624268 98.124132 98.623996 99.123861 99.623725 100.123590
## 197 198 199 200 201 202 203
## 100.623454 101.123318 101.623183 102.123047 102.622911 103.122776 103.622640
## 204 205 206 207 208 209 210
## 104.122505 104.622369 105.122233 105.622098 106.121962 106.621827 107.121691
## 211 212 213 214 215 216 217
## 107.621555 108.121420 108.621284 109.121149 109.621013 110.120877 110.620742
## 218 219 220 221 222 223 224
## 111.120606 111.620470 112.120335 112.620199 113.120064 113.619928 114.119792
## 225 226 227 228 229 230 231
## 114.619657 115.119521 115.619386 116.119250 116.619114 117.118979 117.618843
## 232 233 234 235 236 237 238
## 118.118708 118.618572 119.118436 119.618301 120.118165 120.618029 121.117894
## 239 240 241 242 243 244 245
## 121.617758 122.117623 122.617487 123.117351 123.617216 124.117080 124.616945
## 246 247 248 249 250 251 252
## 125.116809 125.616673 126.116538 126.616402 127.116267 127.616131 128.115995
## 253 254 255 256 257 258 259
## 128.615860 129.115724 129.615588 130.115453 130.615317 131.115182 131.615046
## 260 261 262 263 264 265 266
## 132.114910 132.614775 133.114639 133.614504 134.114368 134.614232 135.114097
## 267 268 269 270 271 272 273
## 135.613961 136.113826 136.613690 137.113554 137.613419 138.113283 138.613147
## 274 275 276 277 278 279 280
## 139.113012 139.612876 140.112741 140.612605 141.112469 141.612334 142.112198
## 281 282 283 284 285 286 287
## 142.612063 143.111927 143.611791 144.111656 144.611520 145.111385 145.611249
## 288 289 290 291 292 293 294
## 146.111113 146.610978 147.110842 147.610706 148.110571 148.610435 149.110300
## 295 296 297 298 299 300 301
## 149.610164 150.110028 150.609893 151.109757 151.609622 152.109486 152.609350
## 302 303 304 305 306 307 308
## 153.109215 153.609079 154.108944 154.608808 155.108672 155.608537 156.108401
## 309 310 311 312 313 314 315
## 156.608265 157.108130 157.607994 158.107859 158.607723 159.107587 159.607452
## 316 317 318 319 320 321 322
## 160.107316 160.607181 161.107045 161.606909 162.106774 162.606638 163.106503
## 323 324 325 326 327 328 329
## 163.606367 164.106231 164.606096 165.105960 165.605824 166.105689 166.605553
## 330 331 332 333 334 335 336
## 167.105418 167.605282 168.105146 168.605011 169.104875 169.604740 170.104604
## 337 338 339 340 341 342 343
## 170.604468 171.104333 171.604197 172.104062 172.603926 173.103790 173.603655
## 344 345 346 347 348 349 350
## 174.103519 174.603383 175.103248 175.603112 176.102977 176.602841 177.102705
## 351 352 353 354 355 356 357
## 177.602570 178.102434 178.602299 179.102163 179.602027 180.101892 180.601756
## 358 359 360 361 362 363 364
## 181.101621 181.601485 182.101349 182.601214 183.101078 183.600942 184.100807
## 365 366 367 368 369 370 371
## 184.600671 185.100536 185.600400 186.100264 186.600129 187.099993 187.599858
## 372 373 374 375 376 377 378
## 188.099722 188.599586 189.099451 189.599315 190.099179 190.599044 191.098908
## 379 380 381 382 383 384 385
## 191.598773 192.098637 192.598501 193.098366 193.598230 194.098095 194.597959
## 386 387 388 389 390 391 392
## 195.097823 195.597688 196.097552 196.597417 197.097281 197.597145 198.097010
## 393 394 395 396 397 398 399
## 198.596874 199.096738 199.596603 200.096467 200.596332 201.096196 201.596060
## 400 401 402 403 404 405 406
## 202.095925 202.595789 203.095654 203.595518 204.095382 204.595247 205.095111
## 407 408 409 410 411 412 413
## 205.594976 206.094840 206.594704 207.094569 207.594433 208.094297 208.594162
## 414 415 416 417 418 419 420
## 209.094026 209.593891 210.093755 210.593619 211.093484 211.593348 212.093213
## 421 422 423 424 425 426 427
## 212.593077 213.092941 213.592806 214.092670 214.592535 215.092399 215.592263
## 428 429 430 431 432 433 434
## 216.092128 216.591992 217.091856 217.591721 218.091585 218.591450 219.091314
## 435 436 437 438 439 440 441
## 219.591178 220.091043 220.590907 221.090772 221.590636 222.090500 222.590365
## 442 443 444 445 446 447 448
## 223.090229 223.590094 224.089958 224.589822 225.089687 225.589551 226.089415
## 449 450 451 452 453 454 455
## 226.589280 227.089144 227.589009 228.088873 228.588737 229.088602 229.588466
## 456 457 458 459 460 461 462
## 230.088331 230.588195 231.088059 231.587924 232.087788 232.587653 233.087517
## 463 464 465 466 467 468 469
## 233.587381 234.087246 234.587110 235.086974 235.586839 236.086703 236.586568
## 470 471 472 473 474 475 476
## 237.086432 237.586296 238.086161 238.586025 239.085890 239.585754 240.085618
## 477 478 479 480 481 482 483
## 240.585483 241.085347 241.585212 242.085076 242.584940 243.084805 243.584669
## 484 485 486 487 488 489 490
## 244.084533 244.584398 245.084262 245.584127 246.083991 246.583855 247.083720
## 491 492 493 494 495 496 497
## 247.583584 248.083449 248.583313 249.083177 249.583042 250.082906 250.582771
## 498 499 500 501 502 503 504
## 251.082635 251.582499 252.082364 252.582228 253.082092 253.581957 254.081821
## 505 506 507 508 509 510 511
## 254.581686 255.081550 255.581414 256.081279 256.581143 257.081008 257.580872
## 512 513 514 515 516 517 518
## 258.080736 258.580601 259.080465 259.580330 260.080194 260.580058 261.079923
## 519 520 521 522 523 524 525
## 261.579787 262.079651 262.579516 263.079380 263.579245 264.079109 264.578973
## 526 527 528 529 530 531 532
## 265.078838 265.578702 266.078567 266.578431 267.078295 267.578160 268.078024
## 533 534 535 536 537 538 539
## 268.577889 269.077753 269.577617 270.077482 270.577346 271.077210 271.577075
## 540 541 542 543 544 545 546
## 272.076939 272.576804 273.076668 273.576532 274.076397 274.576261 275.076126
## 547 548 549 550 551 552 553
## 275.575990 276.075854 276.575719 277.075583 277.575448 278.075312 278.575176
## 554 555 556 557 558 559 560
## 279.075041 279.574905 280.074769 280.574634 281.074498 281.574363 282.074227
## 561 562 563 564 565 566 567
## 282.574091 283.073956 283.573820 284.073685 284.573549 285.073413 285.573278
## 568 569 570 571 572 573 574
## 286.073142 286.573006 287.072871 287.572735 288.072600 288.572464 289.072328
## 575 576 577 578 579 580 581
## 289.572193 290.072057 290.571922 291.071786 291.571650 292.071515 292.571379
## 582 583 584 585 586 587 588
## 293.071244 293.571108 294.070972 294.570837 295.070701 295.570565 296.070430
## 589 590 591 592 593 594 595
## 296.570294 297.070159 297.570023 298.069887 298.569752 299.069616 299.569481
## 596 597 598 599 600 601 602
## 300.069345 300.569209 301.069074 301.568938 302.068803 302.568667 303.068531
## 603 604 605 606 607 608 609
## 303.568396 304.068260 304.568124 305.067989 305.567853 306.067718 306.567582
## 610 611 612 613 614 615 616
## 307.067446 307.567311 308.067175 308.567040 309.066904 309.566768 310.066633
## 617 618 619 620 621 622 623
## 310.566497 311.066362 311.566226 312.066090 312.565955 313.065819 313.565683
## 624 625 626 627 628 629 630
## 314.065548 314.565412 315.065277 315.565141 316.065005 316.564870 317.064734
## 631 632 633 634 635 636 637
## 317.564599 318.064463 318.564327 319.064192 319.564056 320.063921 320.563785
## 638 639 640 641 642 643 644
## 321.063649 321.563514 322.063378 322.563242 323.063107 323.562971 324.062836
## 645 646 647 648 649 650 651
## 324.562700 325.062564 325.562429 326.062293 326.562158 327.062022 327.561886
## 652 653 654 655 656 657 658
## 328.061751 328.561615 329.061480 329.561344 330.061208 330.561073 331.060937
## 659 660 661 662 663 664 665
## 331.560801 332.060666 332.560530 333.060395 333.560259 334.060123 334.559988
## 666 667 668 669 670 671 672
## 335.059852 335.559717 336.059581 336.559445 337.059310 337.559174 338.059039
## 673 674 675 676 677 678 679
## 338.558903 339.058767 339.558632 340.058496 340.558360 341.058225 341.558089
## 680 681 682 683 684 685 686
## 342.057954 342.557818 343.057682 343.557547 344.057411 344.557276 345.057140
## 687 688 689 690 691 692 693
## 345.557004 346.056869 346.556733 347.056598 347.556462 348.056326 348.556191
## 694 695 696 697 698 699 700
## 349.056055 349.555919 350.055784 350.555648 351.055513 351.555377 352.055241
## 701 702 703 704 705 706 707
## 352.555106 353.054970 353.554835 354.054699 354.554563 355.054428 355.554292
## 708 709 710 711 712 713 714
## 356.054157 356.554021 357.053885 357.553750 358.053614 358.553478 359.053343
## 715 716 717 718 719 720 721
## 359.553207 360.053072 360.552936 361.052800 361.552665 362.052529 362.552394
## 722 723 724 725 726 727 728
## 363.052258 363.552122 364.051987 364.551851 365.051716 365.551580 366.051444
## 729 730 731 732 733 734 735
## 366.551309 367.051173 367.551037 368.050902 368.550766 369.050631 369.550495
## 736 737 738 739 740 741 742
## 370.050359 370.550224 371.050088 371.549953 372.049817 372.549681 373.049546
## 743 744 745 746 747 748 749
## 373.549410 374.049275 374.549139 375.049003 375.548868 376.048732 376.548596
## 750 751 752 753 754 755 756
## 377.048461 377.548325 378.048190 378.548054 379.047918 379.547783 380.047647
## 757 758 759 760 761 762 763
## 380.547512 381.047376 381.547240 382.047105 382.546969 383.046833 383.546698
## 764 765 766 767 768 769 770
## 384.046562 384.546427 385.046291 385.546155 386.046020 386.545884 387.045749
## 771 772 773 774 775 776 777
## 387.545613 388.045477 388.545342 389.045206 389.545071 390.044935 390.544799
## 778 779 780 781 782 783 784
## 391.044664 391.544528 392.044392 392.544257 393.044121 393.543986 394.043850
## 785 786 787 788 789 790 791
## 394.543714 395.043579 395.543443 396.043308 396.543172 397.043036 397.542901
## 792 793 794 795 796 797 798
## 398.042765 398.542630 399.042494 399.542358 400.042223 400.542087 401.041951
## 799 800 801 802 803 804 805
## 401.541816 402.041680 402.541545 403.041409 403.541273 404.041138 404.541002
## 806 807 808 809 810 811 812
## 405.040867 405.540731 406.040595 406.540460 407.040324 407.540189 408.040053
## 813 814 815 816 817 818 819
## 408.539917 409.039782 409.539646 410.039510 410.539375 411.039239 411.539104
## 820 821 822 823 824 825 826
## 412.038968 412.538832 413.038697 413.538561 414.038426 414.538290 415.038154
## 827 828 829 830 831 832 833
## 415.538019 416.037883 416.537748 417.037612 417.537476 418.037341 418.537205
## 834 835 836 837 838 839 840
## 419.037069 419.536934 420.036798 420.536663 421.036527 421.536391 422.036256
## 841 842 843 844 845 846 847
## 422.536120 423.035985 423.535849 424.035713 424.535578 425.035442 425.535307
## 848 849 850 851 852 853 854
## 426.035171 426.535035 427.034900 427.534764 428.034628 428.534493 429.034357
## 855 856 857 858 859 860 861
## 429.534222 430.034086 430.533950 431.033815 431.533679 432.033544 432.533408
## 862 863 864 865 866 867 868
## 433.033272 433.533137 434.033001 434.532866 435.032730 435.532594 436.032459
## 869 870 871 872 873 874 875
## 436.532323 437.032187 437.532052 438.031916 438.531781 439.031645 439.531509
## 876 877 878 879 880 881 882
## 440.031374 440.531238 441.031103 441.530967 442.030831 442.530696 443.030560
## 883 884 885 886 887 888 889
## 443.530425 444.030289 444.530153 445.030018 445.529882 446.029746 446.529611
## 890 891 892 893 894 895 896
## 447.029475 447.529340 448.029204 448.529068 449.028933 449.528797 450.028662
## 897 898 899 900 901 902 903
## 450.528526 451.028390 451.528255 452.028119 452.527984 453.027848 453.527712
## 904 905 906 907 908 909 910
## 454.027577 454.527441 455.027305 455.527170 456.027034 456.526899 457.026763
## 911 912 913 914 915 916 917
## 457.526627 458.026492 458.526356 459.026221 459.526085 460.025949 460.525814
## 918 919 920 921 922 923 924
## 461.025678 461.525543 462.025407 462.525271 463.025136 463.525000 464.024864
## 925 926 927 928 929 930 931
## 464.524729 465.024593 465.524458 466.024322 466.524186 467.024051 467.523915
## 932 933 934 935 936 937 938
## 468.023780 468.523644 469.023508 469.523373 470.023237 470.523102 471.022966
## 939 940 941 942 943 944 945
## 471.522830 472.022695 472.522559 473.022423 473.522288 474.022152 474.522017
## 946 947 948 949 950 951 952
## 475.021881 475.521745 476.021610 476.521474 477.021339 477.521203 478.021067
## 953 954 955 956 957 958 959
## 478.520932 479.020796 479.520660 480.020525 480.520389 481.020254 481.520118
## 960 961 962 963 964 965 966
## 482.019982 482.519847 483.019711 483.519576 484.019440 484.519304 485.019169
## 967 968 969 970 971 972 973
## 485.519033 486.018898 486.518762 487.018626 487.518491 488.018355 488.518219
## 974 975 976 977 978 979 980
## 489.018084 489.517948 490.017813 490.517677 491.017541 491.517406 492.017270
## 981 982 983 984 985 986 987
## 492.517135 493.016999 493.516863 494.016728 494.516592 495.016457 495.516321
## 988 989 990 991 992 993 994
## 496.016185 496.516050 497.015914 497.515778 498.015643 498.515507 499.015372
## 995 996 997 998 999 1000
## 499.515236 500.015100 500.514965 501.014829 501.514694 502.014558
Unión en un dataframe
Gráficando los valores
datos_modelo %>%
ggplot( )+
geom_point( aes( x= vector_x, y = E_y_x), col = "blue")+
geom_point( aes( x= vector_x, y = y_observado), col = "red")+
geom_point( aes( x= vector_x, y = y_predicho), col = "brown")+
geom_line(aes( x= vector_x, y = E_y_x), col = "blue")+
geom_line(aes( x= vector_x, y = y_predicho), col = "black")Diferencia entre esperado y el predicho (Error)
## (Intercept) vector_x
## 2.1501694 0.4998644
Construir función
estimar_betas <- function(n,b0,b1,sigma,distribucion){
if (distribucion == "normal") {
errores <- rnorm(40,0,4)
} else
{ errores <- runif(40,-6.928203,6.928203) }
vector_x <- c(1: n)
E_y_x <- b0 + (b1 * vector_x )
#error <- rnorm(n = n, mean = 0, sd = sigma)
#error <- runif(n = n, min=-6.928203, max = 6.928203)
y_observado <- E_y_x + errores
modelo <- lm(y_observado ~ vector_x)
y_predicho <- predict(modelo)
b0_estimado <- modelo$coefficients[1]
b1_estimado <- modelo$coefficients[2]
return(list("b0_estimado" = b0_estimado,
"b1_estimado" = b1_estimado))
}Usando la función
resultados1 <- estimar_betas(40,beta_0,beta_1,sigma,"normal")
# resultados2 <- estimar_betas(40,beta_0,beta_1,sigma)
# resultados3 <- estimar_betas(40,beta_0,beta_1,sigma)
print(resultados1)## $b0_estimado
## (Intercept)
## 3.195765
##
## $b1_estimado
## vector_x
## 0.4823962
Creando vector BETAS
vector_betas <- function(n,beta_0,beta_1,sigma, n_sim, distribucion) {
vector_b0 <- vector()
vector_b1 <- vector()
for (i in 1:n_sim) {
resultados <- estimar_betas(n,beta_0,beta_1,sigma,distribucion)
vector_b0[i] <- resultados$b0_estimado
vector_b1[i] <- resultados$b1_estimado
}
df_resultados <- data.frame(
b0 = vector_b0,
b1 = vector_b1)
return(df_resultados)
}Histograma comparando distribuciones.
histo_b0_normal <- df_resultados1 %>%
ggplot()+
geom_histogram(aes(x = b0))+
ylim(0,110)
histo_b0_uniforme <- df_resultados2 %>%
ggplot()+
geom_histogram(aes(x = b0))+
ylim(0,110)
grid.arrange(histo_b0_normal,
histo_b0_uniforme,
ncol = 2)histo_b0_normal <- df_resultados1 %>%
ggplot()+
geom_histogram(aes(x = b0))+
ylim(0,110)
histo_b1_normal <- df_resultados1 %>%
ggplot()+
geom_histogram(aes(x = b1))+
ylim(0,110)
grid.arrange(histo_b0_normal,
histo_b1_normal,
ncol = 2)